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SOLUTION  OF  TRANSIENT  PROBLEMS 
IN  FREE  SURFACE  HYDRODYNAMICS 

I INTRODUCTION’ 

Lagrangian  methods  offer  the  most  natural  approach  to  transient 
hydro dynamics  problems  which  contain  free  surfaces,  interfaces  or 
sharp  boundaries.  In  practice,  their  use  in  numerical  solutions  has 
generally  been  restricted  to  "well-behaved"  flows,  since  shear, 
fluid  separation  or  even  large  amplitude  motion  produce  severe  grid 
distortions.  The  distortions  arise  from  the  migration  of  mesh  points 
which  were  formerly  neighboring  but  which  have  crossed  or  become 
separated  in  the  flow.  Numerical  solutions  of  the  physical  equations 
differenced  over  such  a mesh  quickly  fail  because  of  the 
inaccuracies  inherent  in  approximations  which  ignore  mesh  point 
rearrangement  and  grid  distortion.  This  problem  is  solved 
by  the  implementation  of  a general-connectivity  grid  in  which 
local  mesh  reconnections  are  made  whenever  the  grid  distorts  appreciably. 

The  techniques  described  in  this  paper  involve  the  use  of  a 
Lagrangian,  finite-difference  mesh  of  connected  triangles  to 
represent  the  fluid  motion,  the  various  interfaces,  and  the  free 
surfaces.  This  approach  is  a fruition  of  efforts  carried  cut  at 
NRL  and  elsewhere  during  the  past  few  years.  Some  of  the  basic 
concepts  were  developed  by  Crowley1  using  the  code  FLAG.  Our  own 

work  has  concentrated  on  the  free-surface  and  physical  consistency 

o_u 

aspects  of  the  problem.  Early  NRL  efforts  centered  on  LIIX'SO  , 
a 2D  MHD  code  with  axial  symmetry.  Attention  more  recently  has 

" .1 

turned  toward  solving  free-surface  problems  in  naval  hydrodynamics. 

Note : Manuscript  submitted  January  17,  1977. 
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A good  example  of  the  type  of  problem  requiring  improved 
numerical  techniques  is  the  problem  of  a large  amplitude  breaking 
wave . '-Jhen  the  wave  top  separates  from  the  wave  due  to  bottom 
shallowing,  nonlinear  wave  interaction,  or  strong  surface  winds,  it 
falls  back  into  the  wave  trough  in  a complicated  flow  which  the 
usual  techniques  cannot  handle.  One  goal  of  this  paper  is  to  present 
Lagrangian  triangular- cell  algorithms  in  which  the  crest  of  the 
numerical  wave  will  be  able  to  separate  from  the  body  of  the  wave, 
fall  under  the  influence  of  gravity  and  any  strong  winds  back 
into  the  trough  of  the  wave,  and  become  reabsorbed  smoothly 
while  following  the  numerical  representation  of  the  fluid  at  the 
surface . 

Figure  1 illustrates  the  use  of  a Lagrangian  triangular  mesh  to 
facilitate  the  solution  of  this  problem.  As  the  neck  of  fluid 
narrows,  the  corresponding  triangle  sides  shrink  to  zero  and  the 
connection  between  the  two  bodies  of  fluid  is  broken.  In  this  way 
the  topology  of  the  problem  has  changed  self -consistently  with  the 
evolution  of  the  system.  Clearly,  any  representation  of  the  problem, 
to  be  adequate,  must  be  capable  of  drastic  local  changes  such  as 
this  to  .reflect  complicated  motions.  It  is  equally  clear  that  a 
numerical  free-surface  capability  which  incorporates  such  reconnec- 
tion algorithms  for  the  solution  of  fully  nonlinear  flews  will  allow 
an  almost  inexhaustible  list  of  important  hydrodynamics  problems 
to  be  solved. 
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Previous  Lagrangian  treatments' ’ ’ have  teen  restricted  to 
simple  flows  or  small  amplitude  motions  because  they  use  a topologi- 
cally rectangular  finite-difference  mesh.  The  major  practical 
problem  with  these  rectangular  methods  arises  from  the  inflexible 

Q 

connectivity  of  the  various  mesh  points.  In  complicated  and 
strongly  sheared  flows,  where  one  element  of  fluid  may  become 
widely  separated  from  a nearby  element,  the  usual  Lagrangian  treat- 
ments break  down  because  a simply  structured  rectangular  mesh 
becomes  too  severely  distorted  to  allow  an  adequate  numerical 
representation  of  the  fluid  flow,  '.'.hen  the  mesh  becomes  sc  distorted 
that  no  greater  Lagrangian  motion  can  be  permitted,  a process  of 
continual  reconing  amounts  to  a form  of  numerical  diffusion.  Thus 
the  usual  Lagrangian  treatments  are  capable  of  extending  linear 
models  significantly  into  the  nonlinear  regime  but  are  not  capable 
of  providing  a satisfactorily  accurate  longtime  representation  of 
highly  complicated  flo\vrs  bordering  on  turbulent  phenomena. 

Another  drawback  of  such  approaches  is  a difficulty  in 
representing  complicated  boundaries  and  structures  because  of  the 
limited  topology  of  the  mesh.  It  is  often  necessary  to  require 
greater  resolution  simply  to  obtain  a satisfactory  initial  grid. 
Rectangular  mesh  approaches  also  appear  to  suffer  a serious  "even- 
odd"  or  computational -mode  instability  which  must  be  overcome  by 
some  form  of  added  numerical  damping.  - ’ This  damping  destroys  the 
reversibility  of  the  algorithm  and  limits  its  usefulness  for  high 
Reynolds-number  flows  even  though  some  care  has  been  riven  to 
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developing  damping  algorithms  which  minimize  this  nonphysical 
diffusion. 

A triangular-element  mesh  has  several  advantages.  The  grid  can 
he  restructured.  Individual  triangles  and  sides  can  he  bisected  or 
rearranged  to  give  new  grid  structures  which  better  represent  changing 
fluid  flows.  Since  the  number  of  triangles  meeting  at  a vertex  is 
variable,  increased  accuracy  in  one  region  of  the  flow  does  not 
force  unnecessary  resolution  in  other  areas  of  the  flow.  This 
versatility  also  permits  both  regular  and  irregular  tessellations 
of  the  x-y  plane  with  triangles.  Triangles,  unlihe  rectangles,  can 
symmetrically  cover  a circle  without  cusps  or  other  local  represen- 
tation irregularities.  Thus  free  surfaces,  complicated  interfaces 
and  boundaries  of  immersed  objects  can  be  represented  accurately  and 
economically. 

The  triangle  is  a much  less  ambiguous  structure  than  a rectangle 
or  higher-order  polygon  and  hence  interpolations  and  integrals  are 
usually  simpler  to  perform.  Moreover,  the  even-odd  problems  of 
rectangular  schemes  appear  to  be  absent,  or  at  least  greatly  subdued, 
in  the  triangular  representation.  Because  of  the  three-sided  figures 
(rather  than  four-sided),  there  is  no  unambiguous  labelling  of  grid 
points  as  even  and  odd.  A vertex,  which  is  one  point  removed  from 
its  neighbor  along  a particular  path,  will  be  two  points  removed  by 
another.  This  dees  not  mean  that  an  even-odd  problem  cannot  occur, 
only  that  it  is  not  a topologically  natural  mode  of  instability  and 
hence  is  considerably  slower  Trowing  than  its  rectan zular-mesh 
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counterpart.  Experience  with  the  SPLISH  code  supports  this;  a 

reversible  algorithm,  or  at  least  a very  weakly  damped  one,  is 

Q 10  1 1 

feasible  using  the  Lagrangian  triangular  mesh  approach.''5  ’ 

Of  course  there  are  also  problems  with  the  triangles.  A 
general  connectivity  triangular  mesh  has  non-trivial  bookkeeping 
problems  associated  with  the  grid  and  its  connections,  the  sides  of 
the  triangles,  and  the  basic  cells  of  the  fluid  dynamic  system. 
Furthermore,  numerical  experience  with  triangles  is  much  more 
limited  than  experience  with  rectangular  meshes.  The  statement  of 
the  hydrodynamic  equations  in  triangular  systems  requires  that 
especially  good  attention  be  paid  to  the  spatial  derivative  terms 
which  are  needed.  There  are,  as  well,  some  intrinsic  numerical 
complications  in  triangular  systems.  The  major  one,  discussed  in 
detail  in  Section  3>  is  the  counting  of  equations  and  free  unknowns. 
This  difficulty  can  be  argued  to  be  the  price  we  pay  for  the  relative 
suppression  of  the  "even-odd"  or  "computational"  modes  observed  in 
more  standard  rectangular  approaches. 

In  this  paper,  we  restrict  ourselves  to  the  study  of  systems  in 
which  the  fluid  is  inviscid  and  incompressible  but  has  variable 
density.  The  conservation  integral  approach  and  definitions  of 
divergence  we  employ,  however,  allow  a natural  extension  tc  compress- 
ible systems.  We  also  restrict  consideration  to  problems  in 

which  the  gravity  is  constant  and  directed  in  the  negative  y 
direction.  These  are  not  necessary  restrictions  but  simplify  the 
analyses  and  allow  the  full  spec trum  of  prcblercs  of  current  interes* 
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to  be  solved.  The  basic  equations  of  the  system  are: 


dV 

p df = * ~P  * psy  ^ C1) 

V . v =0,  (2) 

where  the  fluid  density  o5  pressure  P and  velocity  V,  are  assumed 
to  vary  only  with  x and  y.  Equation  (2),  incompressibility, 
removes  the  sound  waves.  We  will  assume  that  P = constant  or  is 
given  along  free  surfaces. 

The  basic  discussion  of  triangular  grids  and  connection  of  grid 
vertices,  triangle  sides  and  triangle  volumes  is  found  in  Section  2. 
Tessellation  of  the  plane  is  considered  and  the  representation  cf 
complicated  surfaces  and  interfaces  by  a triangular  mesh  is  described. 
Concepts  and  a notation  for  this  triangular  mesh  representation  are 
developed  and  the  vector  aspects  which  can  be  utilised  numerically 
are  given. 

Section  3 is  devoted  to  describing  methods  for  integrating  the 
equations  of  incompressible  hydrodynamics.  A counting  problem 
encountered  in  some  schemes  is  described  and  the  motion  of  the  mesh 
is  discussed.  Two  essentially  reversible  algorithms  are  possible, 
leap-frog  and  centered  implicit,  but  both  require  iteration.  Each 
can  be  coupled  to  any  of  the  several  possible  spatial  treatments. 
Reversibility  is  desired  because  it  reflects  a property  of  the 
inviscid  physical  equations,  and  if  it  is  demanded  of  the  al.g  rithms, 
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major  sources  of  numerical  diffusion  can  be  eliminated.  Throughout 
the  discussion  emphasis  is  placed  on  the  choice  of  algorithms 
actually  implemented  in  the  Cartesian  SPLISH  code.  We  also  discuss 
briefly  questions  of  mass,  momentum  and  energy  conservation  which 
enter  into  the  code . 

Section  k treats  the  important  aspects  of  adjusting  and  rever- 
sibly restructuring  the  mesh.  It  describes  methods  of  setting  up 
the  grid  structure  initially  and  the  various  techniques  which  can 
be  used  to  modify  the  grid  to  better  follow  the  Lagrangian  metier,  o 
the  fluid.  Cr.e  of  the  most  important  aspects  of  this  restructuring 
is  that  strong  shear  flows  can  be  followed  for  long  times  with 
minimal  numerical  diffusion. 

Section  5 is  devoted  to  an  extensive  treatment  of  nonlinear 
free-surfaee  waves  in  which  detailed  comparison  with 
theory  is  made.  Section  6 applies  the  full  numerical  formalism 
with  reconnections  to  the  problem  of  a Kelvin-Kelmholts  unstable 
shear  flow  beneath  a free  surface.  A short  summary  is  contained  in 
Section  7. 

II . Triangular  Mesh  logic 

In  this  section  we  describe  seme  of  the  features  of  triangular 
meshes  and  the  representation  of  the  vertex  interconnections.  The 
techniques  for  adjusting  and  restructuring  the  mesh  during  the 
actual  course  of  a calculation  will  be  attacked  later  in  Section  4. 
Figures  2a  and  2b  show  a section  of  a triangular  mesh  represer.tatic 
with  an  interface  of  a fluid  of  type  I connected  to  a fluid  of  type 
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The  basic  elements  involved  in  the  construction  of  a triangular  mesh 

are  shown . In  figure  2a  a particular  triangle  j , is  shewn  in  heavy 

lines  and  the  various  elements  of  that  triangle  are  labelled.  Three 

'vertices,  V^,  V0,  and  V„,  are  connected  consecutively  by  sides  S1  , 

S9,  and  So . Clearly,  triangle  j shares  these  vertices  and  sides 
— 2) 

with  other  triangles.  The  direction  of  labelling  around  each 
triangle  is  taken  to  be  counter-clockwise  and  the  z direction  is 
out  of  the  page.  Since  the  mesh  can  be  irregularly  connected,  an 
arbitrary  number  of  triangles  can  meet  at  each  vertex.  For  example, 
five  triangles  and  sides  meet  at  vertex  V^.  The  number  of  triangles 
and  sides  meeting  at  a vertex  is  equal  except  near  free-surfaces 
where  there  may  be  no  grid  above  the  surface.  Each  side  is  bounded 
by  two  triangles  (in  general)  and  each  triangle  shares  sides  with 
three  other  triangles. 

Figure  2b  illustrates  several  important  features  of  triangles 
which  are  used  throughout  the  remainder  of  this  paper.  It  is  con- 
venient to  define  a cell  surrounding  a vertex  as  shewn  by  the 
shaded  region  surrounding  vertex  5*  Hie  borders  of  such  vertex- 
centered  cells  are  determined  by  constructing  all  of  the  side 
bisectors  for  each  triangle.  Since  the  three  side  bisectors  all 
intersect  at  a point,  as  shown  for  triangle  j,  there  is  no  ambiguity 
in  constructing  the  vertex  cells  as  indicated.  The  point  of  inter- 
section of  the  side  bisectors,  the  centroid  of  the  triangle,  is  the 
center  of  gravity  for  the  figure  and  this  is  true  in  r-z  as  well  as 
Cartesian  (x-y)  coordinates.  The  three  side  bisectors  of  ar.y 
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triangle  divide  the  triangle  into  six  subtriangles.  These  six 
subtriangles  all  have  the  sane  area  and,  therefore,  each  of  the 
three  vertex  cells  receives  one-third  of  the  area  of  the  triangle, 
'■/hen  a quantity  is  constant  over  a triangle,  the  contributions  of 
that  quantity  from  a triangle  to  each  of  the  three  vertex  cells  are 
also  equal.  These  features  of  side  bisectors  make  the  calculation 
of  cell  volumes  and  cell  areas  particularly  simple. 

Any  computational  representation  of  such  a triangular  mesh 
must  record  all  important  aspects  of  the  mesh  interconnection.  If 
each  vertex,  each  side  and  each  triangle  is  numbered,  lists  of 
interconnections  can  simply  take  the  form  of  an  ordered  series  of 
integers  which  can  be  stored  quite  compactly  in  a computer.  For 
example,  the  information  might  be  stored  in  the  form: 

1.  each  vertex  list  enumerates  (A)  the  vertices  to  which 
the  given  vertex  is  connected;  (b)  the  sides  to  which  it  is 
connected;  and  (C)  the  triangles  which  meet  at  that  vertex; 

2.  for  each  side,  lists  could  record  (a)  the  starting  and  the 
finishing  vertex  of  that  side,  and  (B)  the  triangle  to  the 
right  and  the  left  of  that  side  considered  as  a directed  vector 
from  the  starting  vertex  to  the  finishing  vertex; 

2.  each  triangle  list  could  contain  (a)  the  three  vertices 
of  the  triangle  ordered  in  a counter-clockwise  direction, 

(3)  the  three  sides  cf  the  triangle  between  the  corresponding 
vertices,  again  ordered  in  a counter-clockwise  direction. 
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Since  these  lists  are  not  all  independent,  they  are  not  all 
needed.  Knowing  only  the  vertices  to  which  a given  vertex  is 
connected,  for  example,  still  permits  all  of  the  other  vertex 
information  to  be  determined.  If  vertex  2 is  connected  to  vertex  7, 
a search  of  the  side  lists  can  be  performed  to  determine  which  side 
connects  vertex  7 and  vertex  2, and  thus  the  various  sides  connected 
to  a given  vertex  can  be  determined.  However,  while  maintaining 
redundant  lists  is  costly  in  terms  of  storage,  there  are  real 
advantages  gained  in  increased  computation  speed,  clarity  and  ease 
of  coding. 

It  should  be  clear  that  the  arrays  of  quantities  used  to  define 
the  grid  and  its  motion  involve  the  storage  of  only  local  information 
There  is  no  global  representation  of  the  mesh  in  general,  although 
for  some  specific  geometries,  such  a global  representation  may  be 
possible.  Connection  paths  between  vertices  can  be  determined  only 
by  searching  sequentially  through  neighboring  vertices.  The  lack  of 
a global  representation  in  which  the  vertex  numbering  denotes  a 
corresponding  relative  spatial  position  greatly  complicates  implicit 
calculations.  Poisson's  Equation,  for  example,  must  be  solved  by 
iteration,  and  explicit  time  and  space  derivatives  are  the  best  that 
one  has  any  right  to  hope  for.  Nevertheless,  using  the  incompressibl 
formulation  ensures  that  timestep  limitations  due  to  acoustic  transit 
times  are  not  a problem. 

We  will  use  the  following  notations  and  conventions.  The 
subscript  i will  generally  be  used  to  label  vertices  and  the 
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subscript  j will  generally  be  used  tc  label  triangles.  Thus  II 
denotes  a sun  over  all  vertices  and  I is  a sum  over  triangles. 

j 

The  sum  over  all  three  vertices  of  triangle  j is  denoted  by  II 

-O 

where  the  symbol  is  read  "the  sum  over  vertices  i around  triangle  j" 

Similarly  a sum  over  triangles  around  a central  vertex  c would  be 

denoted  by  S . The  notation  H is  the  sum  over  vertices  i around 
j©  i© 

a central  vertex  c.  In  such  sums  the  sequence  of  vertices  is 
assumed  to  be  counter-clockwise  around  the  central  vertex.  Thus  the 
quantity  A..+^  can  be  used  to  represent  the  triangle  area  of  the 
triangle  with  vertices  (c,  i,  i+l).  Similarly  L.+2_  appearing  in 
H would  be  the  length  of  the  side  radiating  from  vertex  c and 

i© 

separating  triangle  j from  triangle  j+l. 


In  Figure  2 the  area  of  triangle  j is  given  by 


!A.  = (r  - r0)  x (2 


- r 


,) 


(3) 


j ~3  ~2 ' ■'  '~i 

and  thus  h^,  the  height  of  vertex  1 above  the  opposite  side  (S^),  is 
simply 


= (S 


,r0)|/(2A.). 
- J 


00 


Notice  that  the  area  is  a signed  quantity  and  will  be  positive  when 
the  vertices  are  sequenced  in  counter-clockwise  order  around  the 
triangle.  In  fact,  it  is  often  convenient  to  use  the  sign  of  the 
area  to  test  whether  or  not  a triangle  has  inverted  during  the  flow, 
i.e.,  whether  vertex  1,  for  instance,  has  passed  through  side  2. 

As  illustrations  of  these  formulae  and  notations,  the  area  of 
the  cell  centered  on  vertex  i is 


A — aA 

r*. . — — w M.  . 

1 nTTr  .1 
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(6) 


and  the  area-weighted  average  cell  velocity  V,  is  given  by 

v.  = ( Z4a.  v.)/a.  . 

J©  J n- 

Cf  course,  all  the  sums  have  to  be  modified  near  a wall,  an 
interface,  or  a free  surface. 

If  a scalar  function  f is  specified  at  each  vertex  and  is 
assumed  to  be  piecewise  linear  everywhere  between  the  vertices,  the 
vector  gradient  of  f (constant  throughout  the  triangle  and  discon- 
tinuous at  the  triangle  sides)  is  given  by 


(7f)  . 


:1  2A 


-S'  , , 

-L  - 


zx(r2-r. ) 


2A . 
J 


i(2) 


2A . 
J 


For  a vector  field  (V. ] define!  at  all  the  vertices,  the  z 

~i 

component  of  the  curl  and  the  divergence  of  the  vector  are  also 
defined  for  the  triangles  as  follows: 


A.C  = ! (7xV)Z  * dA  = h £ (V.  , + V. ) • (r. 
' J J - ~ ~ J ~ ' i<3)  ~1+1 

j 


A. <?  • V)  . = i (V  • V)  . = dA  = \ (V  +V.  ) x(r  -r.  ) 
J ~ ~ J o ~ ~ j ~1+1  ~1  ~1+1  ~1 


(3) 


fq) 


Equations  (7-°)  define  various  quantities  at  triangle  centroids, 
The  corresponding  quantities  can  usually  be  defined  at  the  vertices 
of  the  finite-difference  mesh  by  area  weighting  the  triangle 
quantities  as  in  Eq.  (6). 
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The  converse  definition  is  also  valid.  A vector  field  may  be 

defined  on  triangles,  where  the  V.  are  constant  throughout  triangles 

‘ '0 

and  discontinuous  at  the  triangle  sides.  For  example,  the  triangle 
velocity  V.  is  given  by 


V. 


(vxV) . = 


V 

L. 


(£i_i"£i+i' 

2A . 

J 


(10) 


where  (ty.  } is  the  vertex-defined  velocity  stream  function  (only  a 
single  z component  in  2D). 

Such  a representation  is  appropriately  both  unique  and  ortho- 
gonal. By  uniqueness  we  merely  mean  that  the  triangle  vector  equation 


V = V0  + V x \lr  (ll) 

~ ~ ~ z* 

uniquely  relates  V and  (0,  ^) (assuming  well-posed  conditions  on 
ar.d  \Ji'  at  the  system  boundaries  where  the  derivatives  are  net 
defined).  Given  i and  on  vertices,  Eqs.  (7)  and  (10 ) allow 
an  unambiguous  determination  of  V at  the  triangle  centroids. 

By  orthogonal  we  mean  that  the  numerical  finite  difference 
operators  satisfy  the  continuum  orthogonality  conditions 

V x V0  = 0,  v • (V  x \|f  ) = 0.  (12) 

The  curl  and  divergence  operators  in  Eqs.  (12)  of  course  act 
on  triangle  quantities  since  t and  v are  assumed  tc  be  vertex 
quantities  as  above  (but  are  otherwise  arbitrary  scalar  fields). 
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The  second  of  Eqs.  (12),  for  example, 

difference  notation  as 

A < V • (V  X \|rj>  = . v • (v 

c 


can  be  written  in  our  finite- 


x i 


(13) 


V. 

1- 


(r.-r  ' 
-1  ~c 


(r. 


i+: 


-r.  ) 


This  is  just  the  flux  of  V x ^ (constant  over  triangles)  out  of 

cell  c.  Clearly  the  coefficient  is  zero  term  by  term  because  of 

c 

the  cross  product.  The  \j/.  and  terms  have  as  coefficients  + \ 

and  - i,  respectively,  because  (r.-r  ) X (r.^-r.)  = 2A.  i.  Thus 

~i  ~c  ~1+1  ~1  1+5 

El-  (13)  is  identically  zero  when  summed  around  an  interior  vertex. 
3y  a virtually  identical  argument  the  first  of  Eqs.  (12)  can  be 
demonstrated  to  be  identically  zero. 

Therefore,  taking  the  divergence  of  V (Eq.  (ll))  expunges  the 
\|x  tern.  Thus 


(r.  . -r. ) 

Vi  ' 44  * 


= V • 70  = 


r 


2, 


zxCr-r.,,)  zx(r.-r  ) zx(r_--r.)- 


(i*0 


+ 0. 


i O’ 


i+4 


i+1  2A.+jk 


c 2A.  ■» 
1+3- 


(r.  -r.) 

~i+l  ~i 


The  second  line  cf  Eq.  (ll) 
operator  and  Eq.  (ll)  is  a 
that  the  coefficient  of  the 
always  negative.  The  natri 


is  the  finite-difference  form  of  the  7^ 
riangular- nrid  Poisson  Equation,  notice 

O 

0 term  is  - H |r.j.-r.  P/lA.+i  and  is 

i©  ~ 

is  diagonally  dominant  and  hence  normal 
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iterative  procedures  for  finding  0 by  inverting  Eq.  (lU)  work  well. 
An  equation  similar  to  (lh)  can  be  constructed,  except  for  sign, 


by  taking  the  curl  of  Eq.  (ll).  Then 

Ac <v  X V>  = V x ( v x '^)  (I?) 

allows  the  determination  of  'i/  from  V. 

III. Finite-Difference  Approaches . 

Given  the  mesh  and  all  its  interconnections,  we  would  like  to 
associate  various  physical  quantities  of  interest  with  the  vertices, 
the  sides  or  the  triangles.  Since  it  is  the  mesh  of  vertices  which 
determines  the  structure  of  the  grid,  it  is  natural  to  assume  that  an 
array  of  positions  is  recorded  for  the  distinct  vertices  in  the 
problem.  To  update  these  positions,  it  is  necessary  to  know  the 
Lagrangian  velocity  of  the  vertices.  Therefore,  we  will  also  assume 
that  the  velocities  of  the  vertices  can  be  found  in  such  a way  that 
integration  of  the  equations  of  motion  for  the  vertex  position  is 
possible  using  these  velocities.  Since  the  velocities  are  assumed 
known  at  the  vertices,  a velocity  field  could  be  constructed  by 
linear  interpolation  within  a triangle  which  is  then  piecewise  linear 
throughout  the  entire  mesh.  This  velocity  field  is  continuous, 
single-valued  and  easily  determined  at  any  point  in  space.  Using 
linear  interpolation,  the  velocities  at  the  center  of  the  sides  are 
simply  the  average  of  the  velocities  at  the  startin'  and  finishing 
vertex  of  that  side. 
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Unfortunately,  the  first  derivatives  of  this  representation  are 
discontinuous  and  the  linear  representation  ensures  that  the  incom- 
pressible equations  of  motion  cannot  be  satisfied  identically  within 
a single  triangle  except  for  a few  simple  cases  such  as  solid  body 
translation.  In  fact,  a real  fluid  element  which  is  initially 
triangular  would  be  very  shortly  transported  in  a real  fluid  to  a 
figure  with  non-straight  sides.  These  deformations  of  the  triangle 
sides  are  a measure  of  the  error  in  the  numerical  methods  we  are 
proposing  and  are  fully  equivalent  to  the  errors  made  in  neglecting 
the  deformation  of  rectangle  sides  in  a topologically  rectangular 
representation. 

The  incompressibility  inherent  in  Eq.  (2)  cannot  be  applied  to 
the  triangles  themselves  for  another,  more  basic,  reason.  The  x 
and  y values  of  the  vertex  positions,  the  variables  available  to 
ensure  conservation  of  triangle  area,  are  fewer  in  number  than  the 
constraints.  This  is  the  previously  mentioned  "counting  problem". 

For  rectangular  grids  a one-to-one  correspondence  can  be  made  between 
quadrilateral  areas  and  the  vertices;  for  example,  associating  a 
quadrilateral  with  its  lower  right  vertex.  For  a triangular  mesh, 
particularly  one  that  allows  reconnections,  no  such  correspondence 
exists. 

In  a quadrilateral  mesh,  the  one-to-one  correspondence  fails 
only  at  the  boundaries.  The  extra  vertex  variables  at  the  boundaries 
are  available  to  satisfy  the  boundary  conditions.  For  a triangular 
mesh 
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or  there  are  less  than  two  triangles  per  vertex.  A one-to-one  corre- 
spondence is  approached  only  for  the  very  restrictive,  and  not  too 
useful,  case  for  which  all  vertices  lie  on  the  boundaries,  when 

'1  = IT  -2.  For  most  cases  of  interest  the  number  of  boundary 
tv 

vertices  is  small  compared  to  the  number  of  interior  vertices, 

and  the  ratio  is  much  worse,  IT,  ~ 2JJ  . (See  the  Appendix  for  details.) 

Although  there  are  always  fewer  than  two  triangles  per  vertex, 
the  number  of  variables  left  free  after  all  of  the  constraints  have 
been  satisfied  are  extremely  few  in  number.  In  this  case,  after  all 
triangle  areas  are  conserved,  the  x and  y positions  available 
could  not  fulfill  boundary  conditions  or  be  used  to  represent  the 
flow  vorticity  with  acceptable  accuracy. 

Once  the  idea  of  conserving  triangle  areas  has  been  abandoned, 
the  next  logical  attempt  is  to  conserve  vertex-centered  cell  areas. 

This  is  attractive  since  the  cell  area  (7  • V)  and  the  cell  rotation 
(7  x T)  are  both  known  for  each  vertex  and  there  are  correspondingly 

two  variables  for  each  vertex,  V and  V , which  are  free  to  satisfy 

x y 

these  constraints.  A great  deal  of  effort  has  been  invested  in 
this  approach  with  mixed  results.  The  cell  area  must  at  all  times 
have  the  value  it  started  with  initially.  The  cell  vorticity 
integral  is  likewise  known.  ’.Then  the  curl  of  Ec.  (l)  is  integrated 
over  a cell  c we  obtain 
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V X V • dA  = 


(17) 


d 

dt 


7P 


X (— ) 


dA 


where  the  term  on  the  right  hand  side  can  be  evaluated  explicitly. 

The  left  hand  side,  the  time  derivative  of  the  average  vorticity 
times  the  cell  area,  can  thus  be  used  to  advance  the  cell  vorticity 
in  time.  Knowing  the  desired  cell  vorticity  and  cell  area,  the  new 
vertex  velocities  can  be  iterated  (ensuring  momentum  conservation 
via  the  equations  of  motion)  until  the  actual  cell  areas  and 
vorticities  are  correct. 

Unfortunately  the  effective  Foisson-lilse  equations  which  result 
are  ncn-local  and  convergence  cf  the  iterations  is  correspondingly 
slew.  While  there  is  no  counting  problem  with  this  approach,  an  allied 
topological  constraint  crops  up  when  mere  than  six  triangles  meet  at 
a vertex.  We  have  been  able  to  show  that  a purely  local  rotation 
conserving  all  cell  areas  is  not  generally  possible  where  seven  or 
more  triangles  meet  at  a vertex  (see  the  Appendix).  In  ether  words, 
more  than  just  the  nearest  neighbor  vertices  must  have  non- sere  velo- 
cities to  conserve  all  cell  areas;  purely  local  vortices  are  not 
generally  possible  and  this  gets  reflected  in  the  convergence  rate  of 
the  iterated  solution.  Similarly,  a purely  vorticity- free  local 
expansion  is  not  generally  possible  when  more  than  six  triangles 
meet  at  a vertex.  This  constraint  is  not  entirely  prohibitive, 
but  the  cost  cf  exact  area  and  vorticity  conservation  with  momentum 
conservation  is  high  in  terms  cf  computer  time  so  faster,  less 
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demanding  algorithms  were  required. 

There  is  no  strict  requirement  to  maintain  physically  conserved 

quantities  any  more  accurately  than  the  truncation  error  levels 

implicit  in  the  finite-difference  algorithms — even  though  such  strict 

conservation  is  esthetically  attractive.  With  this  in  mind  we 

considered  algorithms  which  minimized  a positive  definite  error 

functional  over  the  whole  mesh.  Since  the  desired  curl  {(h)  and 
■* 

divergence  are  known,  and  since  Eqs.  (6),  (8), and  (9)  allow  us 

to  define  at  any  instant  the  actual  values  [C. } and  {Th  },  we  chose  as 
our  error  functional 

S = - c!)2  + W.(D.  - p*)2]  (18) 

The  weights  (uni  and  {W^}  are  free  to  emphasize  small  or  large 
triangles,  fast  or  slow  flow  regions,  etc.  Taking  the  partial 
derivatives  Se/SV.  = 0 gives  as  many  equations  as  one  needs  to 
advance  the  fluid  velocities.  We  never  coded  up  this  last  method, 
however,  because  the  programming  complication  is  prohibitive  and  the 
application  of  boundary  conditions  and  momentum  conservation  is 
obscure. 

The  superficially  attractive  approaches  have  been  mentioned 
briefly  only  to  prevent  (hopefully)  as  many  false  starts  on  the 
part  of  our  readers  as  we  have  undergone.  We  will  now  indicate 
two  distinct  approaches  which  were  successful;  a and  a F-7 
formalism.  Both  define  velocities  as  triangle-based  '.uantities. 
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and  derive  vertex  velocities  from  Eq.  (6).  This  obviates  the 
counting  and  mesh  problems  mentioned  above  through  the  introduction 
of  more  variables. 
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The  first  method  relies  on  the  orthogonality  of  the  curl  and 
gradient  operators  discussed  in  the  previous  section.  We  force  the 
velocity,  now  defined  on  the  triangles,  to  arise  from  a stream 
function  according  to  Eq.  (10).  Such  a flow  is  divergence-free  by 
construction  as  described  earlier;  the  flux  into  or  cut  of  each 
vertex  cell  is  identically  zero.  In  fact,  such  a flow,  as  shown  in 
Figure  4 for  one  vertex  and  its  nearest  neighbors,  is  also  divergence- 
free  on  each  triangle.  Since  V is  constant  throughout  a triangle 
and  since  the  normal  component  of  V is  continuous  across  all 
triangle  sides,  the  flew  field  described  by  V = V x ^ allows  no 
accumulation  of  fluid  in  a triangle  or  loss  of  material  from  a 
triangle.  Figure  3 shows  how  this  can  be  by  displaying  actual 
stream  lines  from  such  a flow.  These  stream  lines  are  closer 
together  where  the  flow  is  faster  and  are  continuous  crossing 
triangle  sides  even  though  their  direction  changes. 

The  introduction  of  triangle  velocities  therefore  makes  it 
possible  to  define  a flow  which  is  divergence-free  on  triangles. 
However,  the  solution  is  not  ideal,  since  the  triangle  velocities 
cannot  be  used  to  advance  vertices,  the  area-weighted  vertex  veloci- 
ties must  be  used.  Therefore  the  triangle  areas  must  necessarily 
change  when  the  vertices  are  advanced,  leaving  only  the  cell  areas 
conserved.  What  remains  is  still  very  attractive  physically, 
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however.  If  the  grid  is  held  fixed,  the  density  cf  marker 
particles  in  the  instantaneous  flow  field  remains  locally  constant, 
while  advancing  the  vertices  leaves  cell  areas  unchanged.  This  is 
the  physically  meaningful  expression  of  incompressibility  which  is 
sustained  where  triangle  velocities  are  the  basic  physical  quantities. 

In  either  formulation  the  incompressibility  assumption  forces 
a physically  nonlocal  character  on  the  flow  and  this  gets  reflected 
in  the  numerical  necessity  of  iterating  a Poisson-like  equation  with 
appropriate  boundary  conditions.  In  the  4-£  form,  the  vorticity  £ 
is  advanced  at  each  cell  during  each  timestep  and  then 
V x (7  x 4)  = £ must  be  solved  for  the  velocity  stream  function  4. 

The  velocity  divergence  is  zero  by  construction.  In  the  primitive 
P-V  formulation,  the  divergence  of  velocity  is  forced  to  zero  by 
choosing  ? to  satisfy  7 • (V  • VV  ) = - V . ^7p,  and  changes  in 
vorticity  are  zero  by  construction.  In  either  case  the  complicated 
grid  connections  virtually  preclude  a fast  direct  solution  or  even  a 
globally  implicit  iteration.  In  solving  these  equations  we  have 
been  using  a simple  point -relaxation  technique  with  a home-grown 
acceleration  formula.  Vie  have  found,  however,  that  simple  extra- 
polation using  the  two  previous  values  of  the  quantity  being 
iterated  (P  or  4')  snakes  the  biggest  impact.  Extrapolation  reduces 
the  residual  by  one  to  two  orders  cf  magnitude  before  iteration  even 
begins.  Vfhat  is  even  mere  important,  however,  is  the  fact  that 
extrapolation  is  most  accurate  for  the  long  wavelength  components 
of  the  solution  desired.  This  is  extremely  helpful  because  the  lon- 
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wavelength  components  of  the  residual  are  the  slowest  to  converge 
under  iteration.  Therefore,  using  extrapolation  not  only  reduces 
the  initial  residual  appreciably,  it  effectively  increases  the  rate 
of  convergence  of  the  iteration  itself. 

The  actual  temporal  integration  scheme  chosen  for  either  formu- 
lation is  reversible,  because  the  physical  equations  being  solved 
are  reversible.'  Reversibility  usually  means  that  artificial 
numerical  damping  is  absent  and  that  higher  Reynold's  number  flows 
can  be  treated  accurately.  Conversely,  the  total  absence  of 
numerical  damping  is  often  accompanied  by  an  even-odd  or  computational 
mode  of  numerical  instability  in  which  decoupling  of  the  finite- 
difference  mesh  variables  occurs  in  some  obnoxious  manner  or  another. 
Previously  we  have  argued  that  our  triangular  grid  algorithms 
minimise  the  virulence  of  this  problem  and  in  Sections  5 and  6 we 
will  show  several  undamped  reversible  calculations.  For  now  we 
concentrate  on  defining  reversible  algorithms  since  damping  can  also 
be  added  where  and  when  needed. 

There  are  two  easy  ways  to  achieve  the  desired  reversibility- - 
using  a leapfrog  x-v  integration  in  time  and  using  a centered 
implicit  _x-v  integration.  Either  approach  can  be  applied  with  either 
the  \jr-£  formulation  cr  the  P-V  formulation  of  the  basic  equations. 

The  choice  appears  to  be  a matter  of  convenience  although  we  have 
used  centered  implicit  in  the  final  versions  of  both  formulations. 

Since  the  flow  is  Lagrangian,  we  can  use  the  following  leap- 
frog time -integration  template 


oo 
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V.  (t+5t/2)  = V.(t-6t/2)  + 6t  a.  (t,  {X},  etc.), 

(19) 

_X.(t+6t)  = X.(t)  + S9  V.  (t+5t/2). 


The  subscript  i labels  the  vertices  (or  the  triangle  centroids. 
The  important  points  are  that  the  acceleration  terms  (a..  } can  be 
calculated  at  time  t using  positions  and  other  quantities  whose 
values  are  already  known  at  time  t and  that  this  "leapfrog" 
integration  is  explicit  as  expressed  above.  Alternately  the 
centered  implicit  formulation  can  be  used. 


V4(t+5t)  =^,'i(t)  + 6t  (t+5t/2), 
X,(t+6t)  = <Xi(t)  + 44  [vjt)  + V.(t+6t)] 


(20) 


where  the  accelerations  _a. (t+5t/2)  are  calculated  either  as  the 
outer  average  i[a,  (t)  + a/t+ot)]  or  the  inner  average  a. (t+6t/2, 
(X^(t)  + /h  ( t+  6t ) } , etc.).  The  formulation  of  Ecs.  (20)  is  clearly 
reversible  and,  insofar  as  the  velocities  V /t+'t)  are  not  known 
a priori,  the  formulation  is  implicit  and  must  be  iterated.  Since 
we  are  iterating  the  acoustic  waves  (incompressibility)  anyway,  this 
ercfcra  iteration  imposes  no  ere at  hardship.  Furthermore,  the  implicit 
nature  sometimes  means  that  longer  timesteps  can  ce  taken  stably. 
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We  will  illustrate  the  \(r-£  method  using  the  easier  algorithm, 
the  leapfrog  method.  We  start  knowing  the  positions  of  the  vertices 
X.  at  time  t and  the  vorticity  at  these  same  vertices  at  time 
t-St/2.  We  also  know  the  parallel  velocity  of  any  free  surface 
vertices  V.  at  time  t-6t/2.  Since  the  pressure  is  constant  at  the 
free  surface,  the  free  surface  boundary  ccndition  can  be  deduced 
from  the  component  of  Eq.  (l)  resolved  parallel  to  the  free  surface. 
If  t denotes  a unit  vector  tangent  to  the  free  surface,  then 


dV 

•f 

"dt 


- gy 


(21) 


Therefore  the  new  free  surface  velocities  have  a parallel  component 
which  is 

V ,(t+6t/2)  = V (t-St/2)  - 6t  gy  ■ r( t).  (22) 


This  is  reversible  and  the  vorticity  can  likewise  be  advanced  tc 
£ (t+6t/2)  reversibly  using  Eq.  (17).  Of  course  the  vorticity  is 
conserved  in  the  Lagrangian  frame  unless  T^  = 0.  .'.ken  To  f 0,  the 
pressure  field  has  tc  be  determined,  a point  we  return  to  shortly. 

The  next  step  is  to  calculate  \!/  from  the  usual  equation 
2,  X (Tx^)  = £.  Here  we  encounter  a snag;  a spatial  grid  at  time 
t+5t/2  is  required  if  \|r.  (t+fit/2)  is  to  be  found.  We  do  this  in 
SPLISH,  cur  '.'’artesian  code,  by  calculating  X.  t+&t/2)  as  the 
average  of  the  new  and  the  old  grid  positions.  Of  course  this  is 
reversible  orovided  we  iterate  exactly  as  in  the  implicit  formula-:  i 
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Cnee  fty.  } has  been  found  the  triangle  velocities  at  time  t-  t/2  can 
be  calculated  using  Z (10  . Forming  the  vertex  area-weighted 
aver  a res  of  V.  from  Zq.  Co)  gives  fv.  ( t+6t/2)  1 which  car.  be  used  to 
advance  the  vertex  positions  X.  from,  t to  t+it. 

The  boundary  conditions  in  she  Poisscr.  Equation  are  single  on  a 
solid  wail?  'i/  = constant.  Cn  the  bottom  of  a fluid  layer  periodic 
in  X,  we  set  \V  = 0.  If  a lid  were  placed  cn  the  top  of  the  fluid, 
the  rc^p^r.t  value  of  y at  the  top  could  be  different  from,  the 
value  cn  the  bottom.  The  difference  in  the  two  constants  would  be 
the  rate  of  flow  across  a vertical  surface  from  the  top  to  the  bottom. 
This  rate  is  independent  of  the  X position  of  the  imaginary  surface 
because  the  fluid  is  incompressible. 

The  boundary  condition  at  the  free  surface  is  merely  that  tire 
component  of  velocity  uarallel  to  the  surface  is  known  for  a free- 
surface  vertex.  The  normal  gradient  of  y,  area-weighted  to  the 
free-surface  vertices,  must  give  the  known  values.  Since  tee  point- 
by- point  iteration  of  the  interior  y points  can  be  performed,  the 
normal  gradient  can  be  adjusted  by  appropriately  choosir.  * the  free- 
surface  -if  values.  An  explicit  evaluation  of  these  values  is 
possible,  Just  as  with  the  interior  points,  but  we  have  chosen  an 
Implicit  formulation  of  this  boundary  condition  which  couples  all 
of  the  surface  \V  values  in  a tridiagonal  system.  In  addition  to 
taking  the  interaction  with  neighboring  surface  vertices  into  account, 
this  method  speeds  convergence  ever  the  whole  mesh,  because 
widely  separated  vertices  can  "communicate  numerical!"”  very  -uickl 
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if  they  are  reasonably  close  to  the  surface. 


r 

When  an  internal  cavity  is  present  in  the  fluid  due  to  the 

existence  of  a bubble  or  to  cavitation,  the  boundary  conditions 

are  greatly  complicated.  For  example,  the  formation  of  a bubble 

has  to  move  the  free  surface  upward  to  conserve  area.  Since 

along  the  surface  is  the  velocity  normal  to  the  surface,  the  net 

d \l/ 

upward  velocity,  ! — d£,  must  be  non-zero.  In  a periodic  system 

OT 

O 

the  formation  of  a bubble  or  cavity  thus  requires  ^ to  be  non- 
periodic even  though  the  derivative  must  be.  In  other  words,  \!/ 
must  have  a branch  cut  whose  jump  varies  with  the  rate  of  expansion 
or  contraction  of  the  bubble.  While  this  nasty  complication  is  not 
insurmountable,  it  certainly  requires  the  calculation  of  a self- 
consistent  pressure  field.  Since  the  fully  nonlinear  source  term 
for  internal  waves,  as  mentioned  earlier,  also  requires  this  pressure, 
a primitive  solution  method  using  P and  V has  definite  advantages 
over  the  \>-'  formulation,  even  though  it  is  intrinsically  a little 
more  complicated  in  2D.  A further  argument  in  favor  of  the  F-V 
formulation  arises  when  we  consider  the  eventual  generalization  to 
three  dimensions.  Then  F is  a true  vector  having  three  components, 
and  three  Poisson  equations  must  be  solved.  In  the  P-Y  formulation 
P would  still  be  a scalar  so  only  one  Poisson  equation  needs  to  be 
solved  even  in  three  dimensions. 


Because  of  these  advantages  the  P-Y  formulation  has  been  the 
major  focus  of  our  work.  All  the  tests  and  examples  ci*-ed  below 
will  use  tills  method,  and  all  further  discussion  will  refer 


exclusively  to  the  primitive  formulation.  However,  it  should  te 
stressed  that  the  formalism  has  also  been  coded  and  tested,  and 
it  has  exhibited  very  good  accuracy  and  convergence.  Although  its 
utility  is  more  restricted,  on  certain  flows  it  is  the  mere  prefer- 
able of  the  two  formalisms  because  of  its  mere  simple  structure  and 
faster  speed  of  execution. 

Our  P-V  formulation  has  P. , V.,  and  X.,  all  soecified  at  times 
~ i ~j  ~i 

7 

t,  t+6t,  t+28t,  ••••.  Using  a split  step  algorithm  , we  integrate 
the  velocities  forward  half  a timestep,  advance  the  grid  one  full 
timestep,  and  then  advance  the  velocities  forward  the  other  half 
timestep. 


V*  = v°  - — (7P)°  - — sr? 

~A  2c,  ■ 'i  2 


(23) 


v;  = i-(v°  + v1}), 

~i  ~i  ~i  ’ 


Xn  = X°  + 6t  V.  , 

~i  ~i  ~i’ 


(2U) 

(25) 


= r ( txn})  • v*, 


V?  = Vf  . Oi.  (7F)u  - — xv 
J 


fit 


\n  fit 


(26) 

(27) 


The  velocity ^ appearing  in  Zq.  (2^)  is  area-weighted  from  'V1?} 
from  Zq.  (2”)  as  calculated  at  the  previous  iteration.  Zqs.  (23' 
and  (2”')  contain  the  evolution  of  (V.l  accordin'  to  t i *•  Laxran'ian 


2 


ecuaticn  of  motion 


Eq.  (26)  is  the  numerical  reflection  of  the  fact 
that  triangle  velocities  must  rotate  and  stretch  as  the  grid  rotates 
and  stretches.  The  transformation  R is  linear  and  given  by  the 
following  three  scalar  equations: 


Tp?  • / v"  ' * { y°  y°  \ 

~2  ” ~1  _ ~j  (^-^l}’ 


# • _ y»)  = V2  • (v°  . X°) 

~J  ~3  ~2  ~3 


(28) 


V*  • (X?  - v")  . 4 • (X?  - X°). 
~i  ~o  ~j  ~l  ~3 


This  transformation  ensures  that  the  vorticity  integral,  calculated 
about  any  interior  vertex,  is  invariant  during  the  displacement  of 
the  grid.  Clearly  Eqs.  (28),  and  hence  Eq.  (2o),  are  fully  reversible. 
It  is  gratifying,  furthermore,  that  the  three  equations  are  not 
linearly  independent  since  only  two  components  of  need  to  ce 
determined. 

As  has  been  established  previously,  the  JR?  and  gravity  terms 
cannot  alter  the  vorticity  either,  since  7 x V?  = C and  .rravity  is  a 
constant.  Cnly  the  physically  correct  variations  of  -±-  car.  cause 
changes  in  vorticity,  exactly  as  they  should.  Thus  the  entire 
algorithm  advances  the  positions  of  the  vertices  and  the  trianvle 
velocities  reversibly  while  evolving  the  correct  vorticity  about 
every  interior  vertex.  The  oressures  {?V}  at  t+it  are  derived  from 
the  condition  that  the  new  velocities  {V?  ^ should  be  diver rence-free 
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at  time  t+6t.  Thus,  from  Eq.  (27)  we  derive  the  pressure  Poisson 


Equation  by  setting 
such  that 


(V 


0.  We  obtain  the  cressures  >F  ^ 

' n ' 


(v  • 


(20) 


The  right  hand  side,  which  can  be  evaluated  at  each  iteration 
explicitly,  is  the  exact  numerical  analogue  of  the  V . (V  • V V) 
term  which  arises  when  the  divergence  of  Eq.  (l)  is  taken. 

Solving  Eq.  (29)  iteratively  for  {F.  } completes  the  calculation 
of  the  timestep.  The  boundary  conditions  on  P for  Eq.  (29)  are 
very  similar  to  those  described  for  the  \|/-£  formulation  earlier.  At 
a free-surface  P = constant  (F  = 0 at  the  top  and  a constant  value 
within  a bubble,  determined  by  the  bubble  Volume,  surface  tension, 
etc.).  A.t  the  bottom  or  on  a wall  cP/5n  must  be  chosen  so  that  the 
velocity  normal  to  the  wall  is  zero.  Thus  an  implicit  tridiagcr.al 
system  can  be  developed  linking  all  pressures  on  the  wall  (or 
bottom)  to  the  nearest  neighbor  interior  values. 

Since  pressures  at  the  boundaries  are  all  specified  through  the 
boundary  conditions,  the  velocities  around  the  partial  cells  at 
these  boundaries  cannot  be  kept  divergence-free  by  iterations  ever 
the  boundary  pressures.  The  boundary  cell  areas  can  be  apportioned 
uniquely  to  the  cell  areas  beneath  the  surfaces,  however,  as  shown 
in  Figure  k.  This  construct  permits  the  pressures  near  boundaries 
to  adjust  for  divergence- free  flow  in  the  expanded  cell  areas, 
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mesh  is  the  flexibility  which  it  permits  to  follow  complicated  flows 
with  a Lagrangian  grid  over  long  times.  To  make  full  use  of  this 
flexibility  requires  that  we  provide  for  several  types  of  adjustment 
and  local  mesh  restructuring  to  maintain  the  uniformity  and  accuracy 
of  the  discrete  mesh  representation.  A mesh  adjustment  is  a non- 
physical movement  or  adjustment  of  the  position  of  one  or  mere 
vertices  which  is  accomplished  without  changing  the  connectivity  of 
the  mesh  vertices.  These  adjustments  are  designed  to  simplify  and 
regularize  the  mesh  and  result  in  the  effective  transfer  of  fluid 
across  triangle  boundaries. 

A restructuring  of  the  mesh,  on  the  other  hand,  dees  not 
generally  involve  movement  of  any  of  the  vertices  but  does  include 
many  sorts  of  triangle  and  side  additions,  reconnections,  and 
subtractions.  In  a sense,  adjustment  and  restructuring  are  ortho- 
gonal procedures--or.e  leaving  the  vertex  positions  unchanged,  the 
other  leaving  the  mesh  connectivity  or  topology  unchanged.  Since 
restructuring  involves  the  changed  position  of  a side,  it  also  can 
involve  the  nonphysical  flow  of  fluid  across  triangle  boundaries. 

Both  adjustment  and  restructuring  represent  departures  from  a 
purely  Lagrangian  description  and  hence  threaten  to  introduce 
unwanted  numerical  diffusion  into  the  system.  To  minimize  the 
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diffusive  ana  ether  errors,  it  is  necessary  tu  pay  strict  attention 
to  mass  and  momentum  conservation  and  to  leave  undisturbed  the 
vertices  and  triangle  sides  lying  along  boundaries,  surfaces,  and 
interfaces.  There  are  a number  of  different  types  of  mesh  adjustment 
and  restructuring  which  are  possible,  too  many  to  cover  in  depth  here. 
Therefore  we  will  fully  develop  here  only  the  most  useful  of  the 
restructuring  techniques,  that  of  reconnection.  We  will  merely  list 
with  a brief  description  the  other  techniques  which  we  have  considered. 
We  tall  conclude  this  section  with  the  illustration  of  a generalized 
mesh  initialization  which  incorporates  many  of  the  adjustment  and 
restructuring  procedures. 

The  triangular  mesh  can  quickly  become  distorted  through  the 
migration  of  formerly  neighboring  vertices  in  the  fluid  flew. 

Without  restructuring,'  the  distorted  mesh  forces  the  computation  of 
derivatives  using  nen-neighbering  vertices,  and  quickly  leads  to 
both  computational  instabilities  and  non-physical  behavior.  This 
situation  is  typified  by  regions  of  long,  narrow  triangles  bordering 
larger  ones.  This  disparity  in  size  also  causes  problems  in  uhat 
timestep  errors  become  severe  because  of  the  short  triangle  sides. 

For  extremely  distorted  triangles,  triangle  inversion  becomes 
increasingly  likely.  Clearly  this  disparity  cannot  be  allowed  to 
increase  indefinitely,  and  some  restructuring  is  needed. 

On  a triangular  grid,  every  nen-boundary  line  uniquely  specifies 
its  two  bordering  triangles.  '.hese  triangles  form  a quadrilateral 
for  which  the  included  line  is  arawn  as  one  of  two  possible  diagonals. 

31 


¥ 


— 1 

Figure  5a  illustrates  a configuration  in  which  the  present  diagonal 
should  obviously  be  reconnected.  In  the  algorithms  employed  her, 
the  shortest  diagonal  for  each  quadrilateral  is  always  chosen,  unless 
reconnection  would  produce  too  large  a disparity  in  triangle  areas. 

This  is  primarily  a safeguard  against  reconnecting  across  inverted 
quadrilaterals,  which  would  produce  negative  triangle  areas,  as  in 
Figure  5b.  It  also  hinders  regions  of  sparse  gridding  due  to  an 
accumulation  of  area  associated  with  a particular  vertex. 

In  order  to  conserve  momentum  locally,  triangle  velocities 
after  a reconnection  must  be  constrained  such  that  the  momentum  of 
the  quadrilateral  is  unchanged.  Furthermore,  to  keep  the  vorticity 
conserved,  choices  of  triangle  velocities  are  further  restricted  to 
those  which  leave  the  vorticity  about  each  vertex  unchanged.  These 
requirements  are  sufficient  to  uniquely  specify  post-reconnection 
velocities  for  the  two  new  triangles.  The  algorithm  resulting  from 
these  constraints  is  reversible.  Replacing  the  reconnected  diagonal 
with  the  original  one  returns  the  triangle  velocities  to  their 
initial  values. 

Further  complications  to  the  algorithms  arise  for  reconnections 
affecting  boundary  vertex  cells  and  by  the  alteration  of  a correc- 
tion term  which  is  carried  to  correct  for  residual  errors  in  the 
pressure  iteration.  Of  course,  general  bookkeeping  of  grid  inter- 
connections must  also  be  updated. 

This  reconnection  procedure  allows  cr.e  row  of  triangles  to 
slip  smoothly  by  another  row  without  having  to  adjust  the  Lagrangian  j 
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vertex  positions  or  velocities.  In  the  corresponding  fixed  grid 
system,  the  triangles  cr  rectangles  bordering  either  side  of  the 
shear  interface  would  soon  become  stretched  unacceptably  and  a 
diffusive  rescue  procedure  would  have  to  be  applied.  The  mesh 
reconnection  described  here  is  a much  less  drastic  change  since  the 
vertex  positions  can  be  left  alone.  The  final  vertex  velocities, 
the  area-weighted  averages  of  triangle  velocities,  wall  be  consistent 
with  conservation  of  vorticity  and  momentum,  and  any  changes  in 
divergence  will  be  resolved  through  future  pressure  changes.  For 
the  special  case  of  an  unperturbed  shear  interface  bordered  by  two 
symmetric  layers  of  triangles,  the  vertex  velocities  would  be  left 
unchanged. 

In  many  cases,  the  Lagrangian  motion  of  the  vertices  naturally 
leads  to  changes  in  vertex  density,  such  as  the  concurrent  accumu- 
lation and  separation  of  vertices  at  a separatrix.  In  such  cases, 
although  the  reconnection  procedures  keep  the  mesh  as  uniform  as 
possible,  the  grid  develops  regions  of  small  and  larve  trial’ vie s , 
and  it  may  become  necessary  to  move  vertices  non-physically.  In 
the  course  of  such  a local  relocation  in  the  fluid  the  total  mass 
and  momentum  of  the  triangles  or  cells  affected  should  of  course 
be  conserved.  Many  possible  formulae  for  the  new  position  have 
teen  tried  ir.  LIJIUS2.  One  of  the  simplest  has  also  proven  to  be 
one  of  the  most  useful--the  new  vertex  position  is  haver,  to  be  the 
average  of  the  surrounding  positions. 


3.3 


In  many  cases  the  non-physical  moving  of  vertices  can  still  be 


avoided  by  adding  and  subtracting  vertices  where  needed.  For 
example,  it  may  be  necessary  to  increase  resolution  in  some  region 
where  the  flow  is  not  naturally  accumulating  vertices.  Such  a 
situation  occurs  in  cyclindrieal  coordinates  when  the  fluid  flow  is 
converging  on  the  axis  somewhere.  Then  the  triangular  zones  become 
larger  and  larger  in  cross-sectional  area,  but  the  radius  of  curva- 
ture becomes  smaller  and  smaller.  Better  resolution  is  clearly 
recuired . 

There  are  at  least  two  ways  of  adding  triangles  to  improve  the 
resolution,  triangle  trisection  and  side  bisection.  In  side  bi- 
section a new  vertex  is  inserted  somewhere  on  a side  which  has 
become  too  long  (for  want  of  another  position  we  can  assume  the 
side  is  bisected).  Lines  from  the  opposing  vertices  to  the  new 
vertex  are  added,  resulting  finally  in  two  new  triangles,  three  new 
sides,  and  a new  vertex. 

Tiie  values  of  physical  variables  at  the  r.ew  vertex  must  be 
determined  by  interpolation,  and  hence  some  numerical  smoothing  is 
implied.  The  improvement  in  resolution  arises  because  the  bisected 
side  becomes  two  sides  and  is  no  longer  required  to  be  straight. 
Thus  increasing  radius  of  curvature  can  be  met  by  increased  reso- 
lution. Furthermore,  there  is  no  restriction  implied  by  material 
interfaces.  If  any  of  the  sides  is  an  interlace,  the  bisection  can 
still  be  performed. 

A new  vertex  can  also  be  added  within  a triangle--not  on  any 
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triangle  side--in  a restructuring  we  call  trisection.  Irisecticr. 
alone  is  irrelevant  because  the  three  new  smaller  triangles  are 
surrounded  by  a single  large  triangle  whose  behavior  is  constrained 
just  as  if  the  new  vertex  were  not  there.  To  be  effective,  tri- 
section  must  be  followed  immediately  by  at  least  a single  side 
reconnection.  The  result  then  is  two  new  triangles,  three  new  sides, 
and  a new  vertex  which  does  not  lie  on  the  original  quadrilateral 
diagonal.  Physical  variables  at  the  new  vertex  are  determined  by 
interpolation  just  as  in  the  case  of  side  bisection. 

Subtracting  vertices  can  be  accomplished  simply  as  the  inverse 
of  these  two  processes.  Per  "inverse  trisecticr.",  a side  is 
reconnected  so  that  the  resulting  configuration  has  a single  triangle 
which  surrounds  a vertex  and  three  subtriangles.  The  interior  vertex 
is  then  erased  and  two  triangles  and  three  sides  disappear. 

Similarly,  a vertex  can  be  relocated  until  two  of  the  sides 
emanating  from  it  form  a straight  line.  If  these  two  sides  belong 
to  triangles  which  share  common  sides,  the  two  interior  sides  which, 
do  not  form  a straight  line  can  be  erased  and  the  two  straight  lines 
combined  into  a single  line.  Here  again  the  starting  configuration 
for  "inverse  side  bisection"  must  be  a quadrilateral  with  four 
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In  s one  cases  the  accumulation  or  separation  of  vertices  may 
be  occurring  far  from  the  region  of  physical  interest.  It  nay 
therefore  become  more  desirable  to  let  the  flew  slip  through  the 
grid  locally.  This  local  Eulerian  behavior  is  similar  to  the  recone 
procedure  employed  by  Chan'  in  a topologically  rectangular  mesh  code. 
Of  course,  one  would  like  to  eliminate  these  numerically  diffusive 
Eulerian  adjustments  within  the  fluid  re  :ion  so  the  numerical  model 
of  the  fluid  mirrors  the  usual  reversibility  characteristic  of  an 
inviscid  fluid  as  nearly  as  possible.  .."her.  Eulerian  slippage  is 
necessary,  however,  the  recently  developed  Flux-Corrected  Transport 

Q 

algorithms'  can  be  adapted  to  minimise  the  diffusion  and  dispersive 
phase  errors  associated  with  the  flow  across  the  cell  boundaries. 

Finally,  there  may  arise  other,  less  general,  mesh  restructur- 
ing which  is  required  by  a particular  problem.  For  example,  periodic 
boundary  conditions  are  often  imposed  to  restrict  a computational 
region.  However,  flow  out  of  the  gcericdic  region  will  carry  t:.e 
grid  with  it.  It  is  necessar:-,  therefore,  to  construct  seme 
means  of  restoring  triangles  which  exit  at  one  boundary  as  triangles 
entering  the  flow  at  the  opposing  boundary.  Since  all  physical 
variables  remain  the  same  for  the  "fictitious’’  triangle  transfer, 
the  procedure  is  mere  bookkeeping,  albeit  complicated,  and  does  net 
alter  the  dynamics  cf  the  flow  in  any  way. 

V.'e  conclude  this  section  with  a brief  description  c;  a general- 
ised mesh  initialization  procedure.  Clearly  ne  soul  sp<  :ify 
hand  all  the  initial  vertex  positions,  all  *d:e  side.',  all  she 
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triangles,  etc.  This  is  both  laborious  ar.d  apt  to  be  repeated  time 
ana  again  as  the  physical  problem  is  changed.  The  process  can  be 
automated  using  general  utility  routines  to  connect  vertices  'and 
store  the  information  in  the  appropriate  mesh  index  lists  , tc 
bisect  sides,  to  reconnect  quadrilateral  diagonals,  and  to  acjust 
interior  vertex  legations. 

This  routine  requires  as  input  an  ordered  set  of  boundary  points 
surrounding  the  desired  region  as  in  Figure  ca,  using  cuts  tc  reach 
an  arbitrary  number  of  interior  subregions  which  are  to  be  excluded. 
The  initialization  program  can  then  proceed  secuentially  around  the 


points,  until  the  region  is  completely  tessellated.  As  shown  in 
Figure  6b,  the  interior  subregions  are  automatically  excluded  through 
the  inclusion  of  the  cuts  in  the  set  of  boundary  points. 

All  quadrilaterals  are  then  repeatedly  scanned,  reconnecting 
diagonals  until  the  shorter  diaaor.al  is  connected  for  each  cruadri- 


longest  remaining  interior  side  is  bisected  repeatedly  until  either 
the  required  number  cf  vertices  has  ceer.  added  cr  all  interior  sides 
are  smaller  than  a specified  length.  Ey  reconnectin'  diagonals 
after  each  bisection,  a fairly  regular  grid  can  be  obtained,  as  in 
Figure  6d.  The  reconnection  and  vertex  adjustment  routines  are 
t:.en  used  iterative! y until  a final  relaxed  mesh  is  obtained  as  ir. 
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This  procedure  can  be  repeated  until  each  separate  region  has 
been  filled  with  an  appropriate  number  of  vertices  and  the  grid 
smoothed.  It  should  be  obvious  that  if  a mere  symmetric  grid  is 
desired  (for  example,  for  use  in  debugging  the  dynamics  of  a program) 
this  same  procedure  can  be  used  to  force  the  grid  to  follow  any 
desired  configuration  with  a minimum  of  effort. 

Since  the  initialisation  begins  with  an  ordered  set  of  boundary 
vertices  and  ends  with  a tessellated  region  corresponding  to  a 
pnysically  identifiable  region,  the  application  of  boundary  conditions 
and  the  storing  of  triangle  and  vertex  physical  variables  can  be 
handled  smoothly.  For  boundaries  which  fall  along  straight  line 
segments,  utility  programs  can  be  employed  sc  that  only  corner 
positions,  vertex  spacing  and  boundary  conditions  need  be  specified 
as  input . 


V . mee-Surface  '..’awes 

Stationary  free-surface  waves  of  varying  amplitudes  in  a 
finite  depth  incompressible  fluid  were  chosen  as  a test  problem  * :r 
the  code  SPLISK.  The  primitive  equation  formulation  was  used,  art 
no  reconnections  were  allowed  to  restrict  the  test  to  the  fur.  ha- 
mental  algorithms.  The  linear  theory  is  well  understc : d ?.r.  rrci'i  h-o 


copious  checks  on  the  numerical  results.  'lonlinear  theory,  th  ■ 
less  developed , also  gives  useful  analytical  results.  'ally, 
ng-time  solutions  are  well  defined  in  the  linear  :ase , thus 


,-v 


yie^cin * a convenient 


of  the  code's 


res 


For  these  tests  a zrid  was  constructed 


represent  a Homo- 


geneous incompressible  fluid  of  finite  depth.  Periodic  boundary 
conditions  were  used  at  the  sides  of  the  region,  permitting  the 
examination  of  an  infinite  wave  train  while  restricting  the  compu- 
tational region  to  one  wavelength.  At  the  rigid  bottom  normal 
velocities  must  be  identically  zero.  In  this  code  this  boundary 
condition  is  met  by  initializing  normal  velocities  to  zero  and 
then  choosing  values  of  F at  the  lower  boundary  which  cause 
normal  accelerations  to  vanish.  The  pressure  at  the  free  surface 
is  taken  to  be  zero. 

The  waves  are  generated  by  specifying  an  impulsive  sinusoidal 
pressure  distribution  at  the  free  surface  at  t = 0.  For  all 
later  times  the  free  surface  pressure  remains  zero.  The  accuracy 
of  the  numerical  solutions  is  tested  through  the  shape,  amplitude 
and  period  of  the  generated  waves. 

A.  Accuracy 

The  linear  theory  yields  an  expression  for  the  period  of  a 

standing  wave  in  a fluid  of  finite  depth  as 

__L 

t = 2n  L gk  tank  ( kh ) „ 2 


where  k - 2tt/\,  g is  gravity  and  h is  the  undisturbed  fluid 
depth.  For  our  tests  \ = 2.5  cm  and  h = 1.  cm  with  g = ~'tC  cm^/sec, 


which  yields  * = 0.12Tl-ii  sec, 


in  the  linear  theory  the  wave  period  is  independent  of  the 
wave  amplitude;  but  the  theory  is  valid  only  wuen  the  rati 


wave  arnlitude  tc  wavelength  and  the  ratio  f -.rave  -unDli  — .de  "■ 


depth  are  bet;,  small 


For  cur  calculations  of  wave  oeriod 


;hese 


ratios  were  A/\  = .0269  and  A/h  = .0672,  since  the  amplitude  had  to 
be  larve  enough  to  visually  determine  the  period.  For  these  values 
we  might  expect  seme  slight  nonlinear  effects  in  the  wave  share. 
However,  as  shown  below,  at  this  value  of  a/\  the  nonlinear  change 
in  period  should  be  negligible. 

Figure  7 presents  the  numerically  determined  wave  periods  for 
different  grid  resolutions.  The  different  mesh  configurations  are 
shown  in  the  inserts  to  the  right  of  the  figure.  To  obtain  a value 
for  the  period  which  had  at  least  four  significant  figures  for  ea^h 
of  these  cases,  without  resorting  to  computer  simulations  lasting 
thousands  of  timesteps,  we  interpolated  between  time steps  for  cur 
time  measurements.  For  a standing  wave,  times  can  be  determined 
most  precisely  as  the  amplitude  passes  through  sere,  when  the  free 
surface  velocities  are  greatest.  'e  defined  the  numerical  value  f 
the  wave  period  by  a least  squares  fit  to  such  time  determinations 
over  computations  lasting  many  wave  periods.  The  error  bars  indica 
the  aggregate  error  in  the  period  due  to  the  uncertainty  in  the 
time  localisation.  It  includes  an  estimate  of  systematic  error  ;.e 
to  a possible  bias  in  selecting  the  precise  time  when  the  amplituue 
was  exactly  nero. 

V.'e  found  the  numerical  period  in  the  limiting  'ase  of  an 
infinitesimally  small  mesh  spacing  by  extrapolation , us in.  • the 
funeti :n 


The  values  of  the  carameters  a and  - are  determined  bv  a least 

c 

squares  fit  to  the  calculated  wave  periods,  r.  The  value  for  the 

oericd  is  ▼ = .12721  = . OOOV  a number  within  '.lc  ner  cent  cf  the 

c 

theoretical  linear  value. 

Contributions  zo  error  in  t from  time step  sine  errors  have 
beer,  minimized  by  reducing  the  timestep  when  the  mesh  size  is 
reduced.  However,  the  timestep  was  not  increased  beyond  ,00U  sec, 
or  roughly  eight  timesteps  per  quarter-period.  In  basically  second- 
order  systems  the  Courar.t  condition  guarantees  that  timestep  error 
terms  are  smaller  than  spatial  derivative  error  terms  because  the 
non-dimensional  timestep  must  be  smaller  than  the  space  step  for 
stability.  This  is  demonstrated  by  the  parabolic  nature  cf  the 
error  curve,  verified  at  six  widely  different  mesh  sizes. 

For  waves  cf  finite  amplitude,  the  wave  profile  should  have 
the  form~^ 


In  the  limit  of  small  ka,  this  expression  reverts  to 
y • h - a ccs  kx,  exactly  t e linear  result.  he  the  reti 
results  then  predict  a smooth  transition  from  a cinu  idal  * : a 


U 


Cur  results 


nearly  trochoidal  form  with  increasing  amplitude. 

also  corroborate  these  predictions,  as  shown  in  Figure  8.  The 

curves  in  this  figure  have  been  normalized  through  the  amplitudes 

of  the  initial  perturbing  oressure  distributions  (A  ).  Inserts  at 

P 

the  right  of  the  figure  illustrate  the  actual  wave  heights  achieved 
for  each  case. 

A more  detailed  comparison  is  shown  in  Figure  9 which  includes 

a superposition  of  the  theoretical  predictions  (Eq.  (31))  and  the 

numerical  results  for  the  A = 8 case  of  Figure  8.  Although  there 

P 

are  slight  differences  in  the  wave  shape,  they  are  net  surprising 
since  Eq.  (31)  is  derived  in  the  limit  of  infinite  depth.  As  the 

inserts  in  Figure  8 show,  the  wave  trough  in  this  case  is  a signi- 

ficant fraction  of  the  depth. 

The  waves  shown  in  Figure  8 were  all  generated  by  an  initial 

sine  wave  pressure  distribution  with  only  the  amplitude  of  the 

impulse  varying  from  case  to  case.  'The  waves  are  labeled  by  the 

impulse  amplitude  (A  ) in  each  case.  The  wave  amplitudes  achieved, 

a , defined  as  ^-(y  - y . ) are  given  in  Table  1. 

n 7 - max  min 

The  linear  theory  predicts  an  amplitude 

a = A^.  -|r-  . (32) 

P P~ 

Using  the  theoretical  value  for  the  period  we  have 

at  = .OI70I  Ap  (33) 

where  we  have  included  am  additional  factor  of  . 2C,  ~c  be  ■v-rplained 
below.  As  can  be  seen  from  the  table,  zhe  numerical!;,  predicted 
amplitudes  deviate  from  the  linear  theory  by  •-.•  percent  at  -..all 
amplitudes  to  7.0  percent  at  lar -e  amplitudes.  ! we-wr.  the  size 


the  mesh  introduces  a finite-difference  truncation  error  in  the 

period,  as  determined  by  Figure  7.  If  we  use  this  value  a for  the 

c 

period,  instead  of  the  value  for  an  infinitely  small  mesh,  we  obtain 

ac  = .Cl68l  Ad  ( 3l) 

as  shown  in  the  table.  The  errors  are  seen  to  be  negligibly  small 

for  the  first  three  cases  and  deviate  from  the  linear  prediction 

only  for  the  obviously  nonlinear  cases.  These  numbers,  coupled 

with  the  observation  that  there  is  no  measurable  increase  in  the 

oericd  up  to  the  A = 8 case  and  an  increase  in  period  for  A = 8 
P p 

and  16.  indicate  that  the  A = 4 case,  used  in  const rcctinc  Fim-re  7, 

P 

is  indeed  a good  approximation  to  the  linear  period.  A further 

check  that  the  percer.ta.re  deviation  for  the  A = 8 and  1'  cases  is 

P 

indeed  due  tc  nonlinear  effects  can  be  made  from  theoretical. 

calculations  of  the  period  for  the  nonlinear  case.  Tadjbakhsh  and 

Keller1^  determine  to  third  order  in  a/\  that  the  period  should 

increase  if  the  initial  fluid  depth  h is  larger  than  C.171.  Ir. 

our  case  the  ratio  h/\  = O.h,  and  sc  the  increase  ir.  period  which. 

12 

is  observed  is  indeed  reasonable.  Penney  ar.d  Friee~“  also  predict 
an  increase  in  their  expansion  to  fifth  order  for  the  case  of 
infinite  depth.  The  change  in  period  predicted  by  each  cf  these 
theories  is  included  in  Table  1,  as  STg,  and  6t  , respectively. 

x Jr 

for  comparison  with  the  amplitude  results.  The  percentage  char.. -re 
should  be  equivalent  since  the  amplitude  achieved  in  ' (27'1  s 

inversely  proportional  to  the  period.  As  car.  ce  seer.  :r  •.  t:.e  ta;  L , 
the  general  trend  of  the  agreement  is  rood,  with,  the  numerical 

13 


results  lying  intermediate  between  the  two  theories 


B.  Stability 

In  numerical  integrations  convergence  to  the  correct  solutions 
is  always  incomplete  and  leaves  a residual  error.  To  ensure  that 
this  residual  error  does  not  accumulate  significantly,  affecting 
both  accuracy  and  stability  cf  the  code,  it  has  been  included 
explicitly  as  "negative  feedback".  For  the  pressure  iteration  in 
which  the  velocity  divergence  around  a mesh  point  is  iterated  to 
zero,  (V^  - V ) / At  is  used  as  a residual  correction,  where  7 is  the 
initial  (t  = 0)  volume  of  the  cell  and  7 is  its  current  volume. 
For  pressures  at  the  bottom  of  the  fluid,  a v source  term  is  added 
to  compensate  for  accelerations  through  pressure  gradient  forces 
perpendicular  to  the  rigid  boundary  which  have  accumulated  due  to 
incomplete  convergence. 


V.’hile  the  feedback  mechanisms  affect  stability  indirectly, 
there  exists  a need  for  more  direct  control  of  numerical  instability. 
As  mentioned  above,  the  divergence  is  sere  only  at  full  timesteps 
and,  therefore,  follows  a time  history  as  shown  in  ligure  10.  .lie 
discontinuous  change  in  divergence  results  from  the  linear  trans- 
formation at  t~2'5t  while  the  linear  change  reflects  the  pressure 
iteration  forcing  the  divergence  to  zero  at  t and  t-ct.  As  seer, 
from  the  figure,  this  scheme  is  marginally  stable  with  the  absolute 
value  f the  divergence  at  e rual  to  that  at  However, 

moving  the  "half"  time step  sli  :htly  forward  cf  center  ensures  that 
the  diver-sense  at  t+  t will  be  smaller  and,  hence , cannot  accumulate 
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secularly,  however,  this  is  dene  at  the  sacrifice  of  perfect  tine 
reversibility  of  the  algorithm.  It  is  the  latter  scheme  which  was 
used  throughout  this  report,  since  many  periods  were  needed  for 
acei;rate  statistics,  for  example,  12  complete  wave  periods  were 
taken  for  the  four  coarser  grids  in  figure  7.  The  forward  shift 
fa  = O.35)  is  also  the  source  of  the  0.37  factor  alluded  to 
previously.  Since  the  impulse  pressure  is  applied  only  for  the 
first  "half"  timestep,  shortening  the  time  also  reduced  the  amplitude 


of  the  impulse. 

As  was  mentioned  above,  long  tine  integrations  are  possible  with 
these  algorithms,  with  runs  typically  500  tinesteps.  for  the  best 
cases,  no  discernible  changes  of  shapes,  amplitudes  or  periods  of  the 
waves  occur  until  t ~ 15"*  Figure  11  illustrates  the  type  of  insta- 
bility presently  found,  here  evident  at  t = 12-.  At  t = 10t  and 
t = 11t  there  is  a slight  flattening  of  the  crests  of  the  waves  pre- 
saging the  mere  violent  onset  of  instability  at  t = 12- . This  cehavio 
reflects  an  acci.mulaticr.  of  truncation  error  which  initially  only 
perturbs  the  grid  near  the  surface,  but  which  grows  in  time  until 
the  grid  is  sufficiently  displaced  locally  to  finally  become  ' our ant - 
unstable.  It  is  not  surprising  that  the  instability  appears,  fat:  er, 
it  is  surprising  that  its  appearance  is  so  retarded  despite  ample 
encouragement  through  the  suppression  of  numerical  averaging  ar.d 
artificial  viscosity.  It  should  also  be  noted  that  there  is  no 
evidence  of  nonlinear  alternating  error  ever,  after  t..e  ir.stabili4"'  ' • 
well  developed,  as  was  su. -rested  for  the  trian.-ular  *r:  :. 

En  summary,  th<  free -surface  wa  e t<  • as  shown  ■ 1 a 
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have  both  excellent  accuracy  and  stability.  The  detailed  examination 
of  its  convergence  has  further  shewn  the  code  to  be  effectively 
second-crder . Finally,  it  has  shown  that  rood  stability  can  be 
achieved  without  artificial  viscosity  or  sratial  numerical  avera.rir.2. 


II.  .he  helyln-'Iehrl'.clts  Instability 

The  dynamics  of  a perturbed  shear  layer  near  a free  surface 
for  an  inviscid,  incompressible  fluid  was  chosen  as  a test  of  the 
numerical  techniques  devised  for  the  reconnections.  The  ilelvir.- 
KeLmholtz  instability  is  an  appropriate  vehicle  for  this  rest  since 
it  has  received  a great  deal  of  attention  both  in  homogeneous  and 
stratified  flows.  It  remains  a fccus  for  research  at  present, 
partially  because  of  the  complexity  involved  in  following  its 
development  even  numerically,  particularly  at  high  Feynolds  numbers, 
and  also  because  areas  su  :h  as  interactions  with  nearby  free 
surfaces  have  remained  almost  totally  unexplored. 

Although  the  purpose  of  this  study  is  to  test  new  numerical 
methods  and  is  net  an  exhaustive  study  of  the  evolution  of  shear 
layer  instability,  a coed  deal  of  insight  int:  the  physical  processes 
involved  can  be  'ained  from  even  coarsely  qridded  simulations  such 
as  these.  In  particular,  the  techniques  presented  here  enable  us 
to  trace  the  evolution  of  the  layer  from  the  initial  perturbations , 
through  their  linear  staves  of  rrewth,  int..  the  manifestly  nonlinear 


growth  of  Kelvin-Keiaholtz  : ill  ws  ai  r ..  into  i te  f nsati  j 

a turbulent  layer. 


Id) 
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A.  Shear  Layer  Dynamics 

The  complete  evolution  of  a shear  layer  near  the  free  surface  of 
a fluid  must  be  calculated  numerically  since  no  analytic,  closed 
form  solutions  exist.  However,  the  early'  growth  should  parallel 
that  predicted  by  linear  theory  until  the  free  surface  becomes 
suff icier.tly  deformed  to  effect  perturbations  in  the  shear  layer . 

The  linear  theory  as  developed  by  Holmbce~~  is  quite  amenable  to 
computer  calculations  and  has  the  further  advantage  of  offering 
additional  physical  insight  into  the  problem.  The  geometry  he 
employs  is  that  of  a finite  shear  lay'er  bounded  by  equal  outer  layers 
of  irrctational  fluid.  He  examines  the  evolution  of  waves  symmetri- 
cally imposed  on  each  of  the  shear  layer  boundaries. 

In  Figure  12  we  have  formed  a symmetric  disturbance  of  the 
layer  by  giving  a small  sinusoidal  perturbation  to  each  boundary' 
of  the  layer.  If  the  boundaries  are  deformed  in  phase,  the  symmetric 
wave  is  in  the  "a-state’’.  For  deformations  which  are  out  of  phase, 
the  wave  is  in  the  "b-state" . Physically,  each  deformation  has 
transferred  vcrticity  across  the  undeformed  layer  boundary,  forming 
a vorticity  wave.  This  wave  in  turn  impresses  a velocity  field  on 
the  surrounding  fluid.  The  subsequent  motion  of  the  layer  is 
derived  from  the  mutual  forcing  of  each  of  the  symmetric  boundary 
waves  through  these  fields,  which  augment  each  other  in  the  b-state 
and  oppose  each  other  in  the  a-state. 


For  perturbations  of  the  shear  layer  with  wavelengths  less 


deformations  travel  with  the  flow  (to  the  right  and  left,  respectively) 
with  oscillating  phase  speed  and  amplitude.  For  wavelengths  greater 
than  1.9  times  the  layer  thickness,  the  deformations  travel  with 
the  flow  from  the  a- state  to  a stationary  state  of  exponential  growth. 
Deformations  reach  the  same  stationary  state  from  the  b- state  by 
propagation  against  the  flow.  Hie  phase  of  the  stationary  state 
is  given  by 


where 
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The  growth  rate,  n,  for  this  state  is 


(^)2  = (l-n  ) (n,  -1) , 
ku  a d 


( 


where  2U  is  the  change  in  velocity  over  the  layer  of  width  i,  and 

k = 2rr/\ . For  \ » d,  kd  « 1 and  n -*  Uk;  i.e.,  the  shear  layer 

behaves  like  a vortex  sheet.  The  wavelength  of  maximum  -rovth 

(4r  = 0)  is  siven  by  \ = 7.9<i. 

dk  ' max 

For  a homogeneous  shear  layer,  the  main  iifferuncvr  • etweer. 
the  behavior  of  a shear  layer  and  a vertex  sheet  is  therefore  a 
short-. ;av'-lor. :ti.  -uteff  for  growth  f tt  initial  perturbation. 


is 


wavelengths  longer  than  the  cutoff,  the  instability  grows  in  a node 
in  which  the  upper  and  lower  boundaries  of  the  layer  remain  phase- 
locked  and  the  growth  proceeds  at  a rate  generally  smaller  than 
that  for  a vortex  sheet.  V.liile  the  vcrticity  associated  with  the 
shear  layer  remains  constant,  its  distribution  changes,  and  becomes 
mere  clustered  as  the  amplitude  of  the  stationary  state  increases 
(Figure  12) . 

If  a statically  stable  density  distribution  exists  within 
the  shear  layer,  however,  spatial  fluctuations  of  particles  within 
the  region  of  density  gradients  will  alter  the  vorticity  baro- 
clinically,  and  the  changing  vorticity  field  -.-Till  produce  a con- 
comitant velocity  field  in  the  surrounding  fluid.  This  field  can 
also  oppose  or  augment  the  symmetric  wave  an  the  boundaries  of  the 
layer.  Yet  the  forcing  is  not  mutual,  since  the  baroclinically 
generated  wave  is  produced  locally  and  is  net  directly  affected  by 
the  deformations  of  the  shear  layer  boundaries. 

The  behavior  of  a stratified  shear  layer  depends  heavily  cn 

the  geometry  of  the  layer  and  the  surrounding  irretaticnal  layers 

and  on  the  functional  form  of  the  density  within  the  layer.  For  a 

1 = 

constant  density  gradient  over  the  layer, Goldstein  found  that  if 
the  outer  layers  were  of  infinite  depth,  the  wavelengths  accessible 
to  the  phase-locked  growth  mode  were  restricted  to  a spectral  band 
which  shifted  to  shorter  wavelengths  with  increasing  F.ic'  ardsen 
number.  Holmcce  found  similar  restrictions  f r t;  is  mode  when 
the  depth  of  the  symmetric  outer  layers  is  h; 
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J)/(J  led) 


> ccth  kh  + tanh 
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(37) 


where  J is  the  Richardson  number, 


Linear  analysis  cf  perturbations  of  shear  layers  in  general 
predicts  a phase-locked  growth  for  unstable  waves,  with  the  growth- 
rates  and  allowed  wave  numbers  determined  by  the  details  of  the 
shear  layer  and  its  stratification.  Although  the  linear  theory  is 
inapplicable  for  later  growth  of  the  layer,  the  principal  mechanism 
for  the  nonlinear  growth  has  already  been  determined. 

As  shown,  in  Figure  13,  at  the  end  of  the  linear  stage  of  growth 
the  vorticity  has  been  concentrated  in  a series  cf  tilted  bands  and 
depleted  from  the  sheared  regions  between  them,  he  will  call  these 
concentrations  of  vorticity  "cores"  and  the  regions  of  strong  shear 
connecting  them  the  "braids”.  A.s  the  vorticity  increases  in  the 
cores,  they  rotate  more  strongly  and  further  stretch  the  braids, 
transferring  even  more  vorticity  to  the  cores.  As  the  roll-up 
continues,  the  portions  of  the  braids  nearest  the  cores  are 
wound  onto  the  cores  also,  accelerating  the  vorticity  transfer. 

For  a homogeneous  fluid,  this  process  could  continue  until  the 
braids  are  totally  depleted  of  vorticity,  leaving  the  core,  in  the 
absence  of  ether  instabilities,  rotating  a' cut  an  a::is  at  the  center 
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For  a stratified  fluid,  the  vorticity  in  the  braids  increases 
with  tilt,  and  stretching  intensifies  the  local  density  gradients 
there.  Vorticity  is  also  generated  in  the  ceres  in  such  a way  as 
to  keep  the  total  vorticity  per  wavelength  constant. 

In  either  case,  the  Kelvin-Helmhcltz  billows  which  are  formed 
increase  in  size  due  to  the  entrainment  of  surrounding  fluid. 
Turbulence  begins  in  the  center  of  the  cere  and  gradually  increases 
as  the  entrainment  progresses.  For  stratified  fluids,  the  braids 
and  the  edges  of  the  cores  become  regions  cf  intense  density 
gradients  along  which  smaller  scale  perturbations  may  appeal'.  Under 
the  combined  influence  of  shear,  strong  density  gradients,  increasing 
billow  size  due  to  entrainment,  and  a growing  central  turbulence, 
the  billows  eventually  coalesce  into  a layer  which  becomes  nearly 
isotropically  turbulent . 

Surprisingly,  at  times  1 ng  in  comparison  with  the  initial 

growth  rate  tim<  s,  ti  is  ti  r'i  lent  layer  is  not  without  structure. 

Micro-layers  appear  w ich,  their  inclination  tc  the  ori  inal 

shear  layer,  car.  be  shewn,  to  be  associated  with  the  original 
If 

billows  , despite  their  subsequent  shearing  and  coalescence. 

Fie  mechanism  for  the  retention  of  the  identity  of  the  original 
billows  long  after  the  formation  of  the  turbulent  layer  has  renaine 
obscure . 


ill 


il  r ii  i i mu  g 
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3.  h'umerical  Solutions 

Figure  ll  illustrates  the  linear  development  of  a shear  layer 
near  a free  surface  as  calculated  cy  SPLISH.  Tlie  inserts  in  th.e 
figure  are  taken  from  the  computer  calculations  and  show  the  location 
of  marker  particles.  The  initial  flow  is  specified  by  the  triangle 
velocities  which  are  chosen  to  be  constant  in  each  layer  of  fluid. 

V.’e  have  given  a velocity  U = 5 cm/sec  (to  the  right)  to  ail  triangles 
above  the  central  discontinuity  and  U = - 5 cm/sec  for  all  triangles 
below  it.  In  terms  of  triangle  velocities,  the  shear  layer  can 
then  be  localized  to  exactly  the  discontinuity.  However,  the  layer 
will  behave  as  if  it  were  one  cell  wide,  simply  because  the  equations 
governing  its  motion  are  differenced  over  a grid. 

A sinusoidal  perturbation  is  given  to  the  central  discontinuity 
at  t = 0 with  an  amplitude  equal  to  2:‘  of  the  fluid  depth  or  20^ 
of  the  shear  layer.  The  perturbation  was  exaggerated  to  help 
visualize  the  flew  in  its  early  stages,  but  because  of  the  rather 
lar re  amplitude, nonlinear  effects  will  be  evident  early  ir.  the 
simulation.  He  effort  has  been  made  to  impose  velocity  fiel  Is 
consistent  with,  these  initial  deformations,  since  th.e  velocities 
will  adjust  themselves  self-consistently  after  a short  period  of 
time.  For  small  perturbations,  these  fludtuations  art  net  lar re, 
and  do  r.ot  significantly  affect  the  dynamics.  It  shouli  he 

emphasised , though,  that  the  initial  c nditi  ns  ar<  i t eigen- 

soluticns  for  a riven  layer  pertv.r’;  at  ion,  '•*  morel'  re  ore : m* 
an  instantaneous  deformation  at  tie  -outer  f • t aver. 


a symmetric  wave 


As  shc/n  in  oirure  lb. 


near  tr.e  "-state  is 


imoressed  cn  the  shear  layer  by  tr.e  central  perturbation.  Initia 
the  symmetric  wave  is  small  compared  to  the  central  deformation, 
phase  speed  of  the  symmetric  wave  is  reversed  by  the  velocity  fie 
associated  with  the  central  wave  as  it  passes  through  tr.e  b-state 
toward  the  a-state.  The  symmetric  wave  then  enters  the  predicted 
phase-locked  mode  of  growth  as  it  nears  the  stationary  state.  As 
shown  in  the  graph,  the  numerical  calculation  determines  the 
staticnarv  state  to  be  at  m = - 7C.C-  and  wives  a rrovth  rate  of 


s 

n = 19.6.  From  Iqs.  (35)  and  (36)  for  d = 0.1  cm,  X.  = 1.0  cm  and 

o 

U = 5 cm/ sec,  the  linear  theory  predicts  cy  = - b7.5  2nd  n = 1". 
in  quite  good  agreement  with  the  numerical  result.  Roughly  hal 
of  the  small  discrepancy  is  cue  to  nonlinear  effects,  evident  by 
the  peaking  in  tr.e  last  insert,  and  the  rest  is  due  to  the  finite 
resolution  of  tr.e  grid. 

To  check  whether  this  agreement  was  merely  fortuitous,  the 
simulation  was  rerun  with  half  the  initial  perturbation . Tr.e 
numerical  results  were  c = - • . and  r.  = 1 .2.  even  closer  t 


s 

the  predictions  of  the  linear  theory,  as  should  be  errpected. 

We  have  also  performed  this  test  ir.  the  limiting  case  c: 
letting  the  perturbation  equal  sere.  Ir.  this  case  the  shear  lave 
is  in  an  unstable  equilibrium,  and  any  error  in  the  assignment  o: 
post-re connect!  n "el  :ities  should  be  mirrored  in  a similar] 
exponentially  gr  wing  listur  ance,  lis  I a:  een  run  : r 

tinesteos  witl  ihanges  in  the  n*i  as!  fron  1 Lnil  al 
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shear  (vertical  grid  positions  indicate  errors  £ 10  ' fcr  a single 


precision  calculation) . The  final  grid  looks  exactly 


y as  iz  aid 


initially,  yet  each  vertex  in  the  irrctaticr.al  fluid  above  the  shear 
layer  has  traversed  the  entire  grid  ten  times. 

Figure  If  presents  another  check  on  the  accuracy  of  the 
computer  calculations  in  the  linear  regime . Here  the  wavelength 
of  the  computational  region  is  increased  to  \ = 1.5  cm  instead  of 
\ = 1.0  cm  as  in  the  previous  simulations.  Although  the  initial 
perturbation  has  this  same  wavelength,  the  half -wavelength  wave, 
with  X = 0.75  cm,  has  a substantially  larger  growth  rate.  This  is 
because  X = 7.5d  and,  as  shown  above,  this  is  quite  near  the  wave- 
length for  maximum  growth  fcr  the  layer,  X = 7.Qd.  Physically,  we 
would  expect  that  such  a layer  would  develop  through  the  half- 
wavelength  mode  due  to  small  random  perturbations.  Numerically,  we 
expect  the  same  behavior,  since  the  initial  perturbation  is  non  an 
eigensolution  for  the  layer. 

As  shown  in  Figure  15,  the  symmetric  wave  that  initially 
develops  is  in  the  b-state  (t  = 0.02),  but  here  it  continues  to  -row- 
in  the  b-state  at  the  expense  of  the  central  wave.  r-y  t = 0*0?. 
the  central  perturbation  is  reduced  to  nearly  one-tenth  of  its 
original  amplitude  and  has  departed  from  its  sinusoidal  form, 


whereas  the  symmetric  wave  is  well  developed  and  is  predominantly 


in  a full-wave length  (\  = 1.5)  b-state.  As  tk 


ne  symmetric  wave 


moves  further  off  the  b-state,  the  fastest  trowing  per tu rbat i or.s 


of  the  central  wave  are  not  the  full-traveler- " modes  but  the  - alt'- 


•.■rave length  (t  = 3.12).  By  t = 0.16,  the  growing  central  pertur- 
b at  ion  has  induced  a half-wavelength  node  cn  the  symmetric  wave, 
and  at  t ~ 0.22  the  layer  is  clearly  predominantly  perturbed  by  a 
half-wavelength  disturbance.  The  symmetric  wave  at  this  time  is 
net  far  from  the  stationary  state  for  the  half-wave length  mode, 
and  its  rapid  growth  clearly  carries  the  layer's  development  into 
the  nonlinear  regime  in  just  0.08  sec  more  (t  = 0.30). 

Figure  16  presents  the  results  of  calculations  for  the  nonlinear 
growth  of  a shear  layer  into  a Kelvin-Helmholtz  billow.  This 
simulation  is  a continuation  of  the  linear  stages  shown  in  Figure  ll. 
The  lower  four  inserts  present  vorticity  contours  at  the  same  times 
as  the  marker  particles  displayed  above.  Since  the  initial  shear 
layer  for  this  problem  is  one  cell  wide  with  constant  velocities 
above  and  below  the  layer,  those  cells  defining  the  layer  will  have  a 
uniform  vorticity.  In  the  inserts  the  vort.icity  is  contoured  over 
the  mesh  triangles.  The  total  width  of  the  layer  therefore  arrears 
twice  as  large  since  the  contours  extend  to  neighboring  vertices. 
Because  the  algorithms  identically  conserve  vorticity,  those  cells 
initially  in  the  layer  always  have  the  same  vorticity,  but  the 
distribution  of  these  cells  will  change. 

In  the  first  insert  the  vorticity  has  begun  to  cluster  into 
bands  with  an  accompanying  thinninv  and  stretching  of  vorticity 
outside  the  hands.  Due  to  cur  finite  resolution,  the  thinning  : ■ 
this  plot  is  evidenced  more  by  the  iensity  of  vcrtici^y-  •arrovn-- 
- gy  • gy  acisua  ■ r ■ • : ■ f ■ • . • . 
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enlarged  band  tilts  and  breaks,  entraining  the  surrounding  fluid. 

The  roll-up  and  gathering- in  of  vorticity  is  shown  by  the  core 
rotation  as  well  as  by  the  movement  of  two  vorticity- carrying  cells 
left  in  the  braids.  In  the  last  frame  the  Kelvin-Helmholtz  billow 
has  matured  and  continues  to  wind  in  the  remaining  vorticity. 

The  marker  particles  carry  complementary  information.  The 
nonlinear  peaking  at  the  end  of  linear  growth  rapidly  leads  to 
breaking  as  the  band  tilts.  In  finer  resolution  this  breaking  at  the 
band  edges  would  be  represented  by  small  vertices.  In  simulations 
using  smaller  initial  perturbations  these  vertices  are  replaced 
by  a much  more  uniform  roll-up  about  the  central  core.  Intrainment 
in  those  cases  is  always  directed  toward  the  center  of  the  core. 

The  smaller  scale  motion  shown  here  is  mere  effective  in  entraining 
fluid.  The  billow  grows  mere  rapidly  through  the  increased 
transfer  of  fluid  to  the  core  by  the  additional  entrainment  of  the 
small  vertices. 


In  these  clots  the  stretching  and  thinning 


tne  orards  i: 


much  mere  evident,  particularly  as  the  billow  matures  and  continues 
to  wind  the  braid  around  the  cere.  Tice  oinks  in  the  braids  at  the 


locations  of  the  remaining  vorticity  are  no 


t pr.ysical,  cut  are  . i.e 


finite  resolution.  However,  they  do  not  represent  a numerical 
instability,  since  these  perturbations  do  net  increase  with  tire. 

The  details  of  these  simulations,  particular!;,  '.hose  wi* .. 
smaller  initial  perturbations « are  reminiscent  of  the  results 


obtained  by  Fatnaik  et  al.  in  their  1 autifull  res 


for  which  at  higher  Reynolds  numbers  the  cere  rotation  becomes  rapid. 
The  primary  differences  appear  to  be  the  result  of  our  nearly  four 

times  coarser  'ridding.  For  example,  ncr.e  of  their  simulations  has 

« 

a well  developed  craid  deformation,  even  at  their  highest  Reynolds 
numbers,  which  confirms  that  their  origin  here  is  due  to  finite 
resolution. 

Of  course  the  main  difference  between  these  calculations  and 
previous  ones  is  the  presence  of  a free  surface  near  the  shear  layer. 
In  Figure  17,  the  evolution  of  the  free  surface  near  a h-H  billow 
is  presented  as  a function  of  time  for  two  different  amplitudes  of 
the  initial  perturbation.  lime  begins  at  roughly  the  end  of  the 
linear  growth  stage  and  ends  at  the  formation  of  a mature  billow 
in  each  case.  The  difference  in  the  times  is  accounted  for  by  the 
larger  initial  amplitude  in  the  second  case  and  by  the  more  rapid 
development  of  its  billow  due  to  the  small-scale  motion. 

The  general  shape  of  the  free  surface  is  the  same  in  both 
cases.  The  free  surface  is  lowest  iirectly  a;  • vc  the  billow  and 
highest  ever  the  braids.  Wave  steepening  Ls  Lriven  by  th«  : 1 id 
entrainment  into  the  billow  below.  he  steepest  portion  the 
surface  wave  remains  above  the  strongest  vertical  velocities  of 
the  entraining  fluid,  movin'  to  the  riot  as  the  '.  ill  w -rows . 


the 


Ap  = 0.1  case.  This  occurs  because  the  billow  builds  more 
rapidly  for  the  larger  perturbation,  and  this  mere  rapid  growth 
brings  the  billow  closer  to  the  free  surface. 

Throughout  this  section  ail  examples  have  used  homogeneous 
shear  layers,  and  connotations  such  as  "braids"  have  teen  used  cn 
the  basis  of  path  lines  for  Lagrangian  particles  alone,  hith  the 
addition  of  stratification,  such  entities  become  loci  for  changing 
density  gradients,  and  may  in  themselves  now  baroclinically  generate 
vorticity.  Figure  IS  presents  vorticity  contours  for  a stratified 
shear  layer  using  an  even  coarser  resolution.  At  t = 0.1c,  vorticity 
generation  in  the  braids  is  already  evident,  as  well  as  the  more 
obvious  vorticity  depletion  in  the  regie n which  will  scon  develop 
into  a core.  As  the  layer  rolls  up  at  t = C.2h,  this  trend  is  quite 
striking,  and  the  small  extent  of  the  core  is  emphasized.  "• 
t = 0.2v,  the  core's  vorticity  is  severely  depleted  and  the  core 
is  winding  in  the  large  concentrations  of  vorticity  in  the  braids. 

At  t = 0.32,  the  regeneration  of  vorticity  in  the  thinned  l rails 
is  already  well  underway. 

It  is  clear  that  despite  the  addition  of  density  stratification, 
the  shear  layer  is  still  unstable  to  phase-locked  rcll-up.  This  is 
exactly  what  is  predicted  by  So.  (37);  the  limits  are  l.h  • and  3.  •, 

and  the  value  of  (l  + j)./(kd/2)  is  2.06. 

he  checked  to  see  if  stabilization  could  be  achieved  with  a 
constant  gradient  by  increasing  the  depth  of  -:.e  s.  ear  layer  v 
l * ^ and  Ap  to  3,2.  lere  Eq.  37)  predicts  a strongly  stabiliz  J 
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layer 


Ihe  numerical  results  conformed  to 


• for  \ = 1.3  disturbances . 
this  expectation.  In  the  previous  case  the  growth  rate  for  the 
perturbation  cad  merely  teen  attenuated  by  the  stratification.  Here 
ti.e  deformation  begins  collapsing  immediately,  and  at  t = 0.-  is 
nearly  gone.  The  layer  continues  its  progress  through  the  mean 
position  and  undergoes  a standing  oscillation. 

The  computer  simulations  therefore  seem  to  agree  in  all  respects 
with  predictions  of  both  linear  and  nonlinear  theory.  However, 
there  is  the  added  advantage  in  that  the  simulations  can  be 
extended  economically  to  longer  times.  In  fact,  we  have  performed 
calculations  at  the  coarser  resolutions  which  have  followed  the 
shearing  of  billows  into  their  coalescence  with  neighboring  billows. 

A sample  calculation  is  illustrated  in  Figure  19.  .At  t = 0.35'? 
the  entire  outer  edge  of  the  billow  is  ringed  by  former  braid  wound 
onto  the  billow,  and  in  a stratified  fluid  would  be  a locus  for 
high  density  gradients.  Frier  to  this  time,  the  billow  has 
already  sheared  and  irrotatienal  fluid  has  beer,  trapped  between 
the  spreading  billow  and  the  braid,  accentuating  the  stretching  of 
the  braid  even  mere.  As  shewn  at  t = 0.29c  and  t = O.-lf,  the 

sheared  billow  continues  to  overlap  the  braid  in  this  region,  and  by 
t = 0.^3'c,  yet  another  layer  of  irrotational  fluid  is  about  to  be 
pressed  avainst  this  already  twice  overlaid  portion  of  the  braid. 
Although  the  marker  density  is  suite  small  at  later  times,  it  is 
possible  to  see  that  at  t = and  t = this  portioi 

b rai  i is  still  intact  ar. : further  stretched,  being  mere  1;  verlai 
by  the  shearing  billow  In  successive  layers  cf  turbulent  and 
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irrotational  fluid. 


This  development  is  mirrored  in  the  vorticity  contours.  It 
is  immediately  obvious  that  the  vorticity- carrying  vertices  perform 
localised  rotations  during  the  shear,  and  do  net  migrate  horizontally. 
Ey  transposing  an  image  of  the  braid  locus  onto  the  contours  and 
tracing  the  movement  of  individual  vertices,  it  can  be  seen  that  the 
braid  effectively  divides  the  fluid  into  separate  regions,  despite 
the  increasingly  isotropic  turbulence  throughout  other  regions  of 
the  sheared  billow. 

The  long-time  simulations  therefore  indicate  that  portions  of 
the  braid  survive  well  into  the  turbulent  regime,  despite  the 
fragmentation  of  the  billows  which  created  them.  These  extremely 
thinned  and  overlaid  braids  inhibit  motion  across  themselves,  and 

serve  as  an  excellent  locus  for  the  building  of  the  micro-layers 

• 16 

seen  xn  late-time  turbulent  regions  as  noted  by  Thorpe.' 

VII . Conclusions 

In  this  paper  we  have  described  an  entire  series  of  techniques 
and  algorithms  for  using  a Lagrangian  £D  mesh  of  triangles  to 
represent  and  solve  free  surface  problems  in  incompressible  hydro- 
dynamics. The  results  of  our  research  in  this  area  have  been 
incorporated  into  the  free- surface  hydrodynamics  code  SPLISH, 
which  was  used  in  calculations  to  test  the  accuracy  and  stability 
cf  the  algorithms.  These  results  on  the  problem  of  free-surface  waves 
have  shown  that  the  code  is  an  extremely  accurate,  second-crier  cede 
and  has  good  convergence.  It  also  has  the  ability  to  calculate 


tiO 
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viscosity  or  numerical  averaging. 

sir^ls/ticns  o i'  "tr.G  '"0 2_vi. n— " '0 j-CC1  cZ.cz  -C3.V0  s ■ c~.m 

that  it  is  indeed  possible  tc  devise  accurate  mesh  reccnr.eaoicr. 
algcrithms  ■which  y'ermit  long-time  solutions  even  in  the  presence 
of  strong  shear.  Here  too,  the  numerical  solutions  have  teen  found 
to  be  accurate,  agreeing  with  both  linear  and  nonlinear  theory  for 
both  hone aeneous  and  stratified  flew.  It  has  been  encouraging  that 
even  coarsely  rridded  problems  such  as  these  not  only  preserve  the 
details  of  previous  higher  resolution,  lulerian  cal  relations , but 
even  elucidate  current  problems  such  as  the  behavior  of  free 
surfaces  near  a shear  layer  and  the  development  of  microlayers  in 
turbulent  regions  formed  by  coalescing  :ie Ivin- Helmholtz  billows. 

Our  experience  with  this  code  is  cv  no  means  extensive,  however 
Prior  development  through  !HX'S2,  a triangular  code  for  electro- 


such  as  TCRUS2,  an  r-c  quasi-static  Mil  cede  for  calculating 
plasma  equilibria  in  Tokamaks , have  led  to  further  experience  in  hew 
tc  difference  properly  equations  ever  a triangular  mesh.  Yet  in 
all  cases  we  have  preceded  only  with  those  schemes  which  seemed  mcs~ 
productive  fer  the  immediate  future.  In  this  paper  we  have  tried 
to  illuminate  seme  of  the  pitfalls  which  we  did  net  avoid  in  the 
hope  that  others  will,  it  is  only  through  cumulative  experience 
such  as  this  that  an  understanding  cf  '■.•.•arc  iynamics  usin  * Lagranviar. 


:riar. vuiar  meshes  will  begin  tc  be  ccmmiete. 
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Aopendix 


Tct30 logical  Constraints 


Since  little  is  know,  about  the  convergence  of  lagrangian 
algorithms  differenced  over  triangular  meshes,  it  is  important  tc 
obtain  a feel  for  the  restricticr.s  that  the  grid  itselr  imposes  or. 
the  physical  solutions  tc  a problem.  Perhaps  the  most  fundamental 
restriction  is  the  relation  between  the  number  cf  vertices  and  the 
number  cf  triangles,  the  "counting  problem"  mentioned  in  the  tent 
and  derived  in  detail  below.  Other  constraints  may  be  more  subtle 
such  as  the  case  presented  here  in  which  local  grid  structures  can 
force  non-local  behavior  onto  the  physical  e :uaticns  which  are 
differenced  over  the  mesh. 

Ic  define  the  "counting  problem"  we  '■  n ::r.siier  first  the  'as 
ex'  9,  cissn. nTuILs-i’  3.'tz.on&2_  ^0  ri  ri  ?_ r s , ..t.  1 . i *,  9..  Tr.c 


icier.  :e 


derivation  was  pointed  out  to  us  by  lr. 

Applications,  Inc.  he  want  tc  find  the  sum  of  Les  sut ten 

by  each  vertex.  Per  an  interior  vertex,  we  have  simply 


< 

0 


. or  vertices  cr.  tr.e  counaary,  out  net  at  a vrner 


V .v  _ _ 
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winle 


, icr  corner  vertices, 


Z k = 


rr/? 


Therefore,  the  total  is  the  sun  of  Eqs.(A-l) 


( A 


1 $ = (N  -II.  )2tt  + (K  -U)n  + Mtt/2), 


V V,  ' ' V, 

b o 


where  IT  is  the  total  number  of  vertices,  and  '*  the  nur.be ^ c 
v vb 

boundary  vertices.  Simplifying 


* = 2n 


V 


This  sum  can  be  calculated  another  way,  through  the  triangles  : 
the  mesh, 


where  IT  . is  the  total  number  of  trianrles.  cmuatine  T's. 

O w - ' 

(A—:)  gives  the  result, 
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The  equations  in 
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able  I.  Comparison  oi  numerical  results  with  linear  ar.d  nonlinear 
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Fig.  4 — Enlarged  cell  areas  near  a boundary.  Because  pressures  at 
boundaries  are  fixed  by  boundary’  conditions,  the  pressures  of  ver- 
tices adjacent  to  the  boundary  are  used  to  keep  the  additional  sur- 
face areas  divergence-free  as  well. 
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RECONNECTION  OF  QUADRILATERAL  DIAGONAL 


DIAGONALS  FOR  INVERTED  QUADRILATERAL 


Fig.  5 — Portions  of  a grid  illustrating  possible  reconnections.  In  part  (a) 
the  dashed  diagonal  will  be  chosen  for  the  shaded  quadrilateral  rather  than 
the  present,  longer,  diagonal.  In  part  ( b>  the  diagonal  cannot  be  recon- 
nected since  the  alternative  diagonal,  though  shorter,  lies  outside  the  “in- 
verted" quadrilateral. 


Richardson  extrapolation  of  the  wave  period  as  a function  of  the  grid  spacing. 
Inserts  illustrate  the  variation  in  mesh  sizes  used. 
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Fig.  8 — Transition  of  the  wave  profile  from  the  linear  to  the  nonlinear  regime. 
The  wave  profiles  are  normalized  by  the  amplitudes  of  the  initial  pressure  per- 
turbations P = Ap  sin  kx.  Note  the  compressed  scale  in  the  x-direction  for  the 
inserts. 
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hie.  !)  Comparison  of  tlx*  numerical  results  with  the  fifth-order 
theoretical  prediction  in  the  limit  of  infinite  depth 


o 


TIME 


B)  FORWARD  TIME  SPLITTING 

Fig.  10  — The  time  development  of  the  divergence 
for  two  differencing  schemes 
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END  OF  LINEAR  GROWTH 


WAVE  BREAKING 


ROLL-UP  AND  ENTRAINMENT  BILLOWS  WITH  CORE  TURBULENCE 

Fig.  13  — The  formation  of  a Kelvin-Helmholtz  billow  from  a perturbed  shear 
layer.  In  the  homogeneous  case  the  vorticity  in  the  shaded  area  remains  constant 
throughout  the  different  stages  of  growth.  For  stably  stratified  layers,  vorticity 
is  generated  and  depleted  baroclinically  along  the  deformed  layer  according  to 
its  tilt. 


Fu;.  la  — Growth  of  a favored  harmonic  of  the  initial  perturbation 
Here  cl  = 0.1.  A0  = 0.01.  6t  = 0.001,  and  the  mesh  is  11  X 11.  The 
initial  perturbation  of  wavelength  X = 1.5  iias  a theoretically  predicted 
growth  rate  of  n = 15.4.  Its  first  harmonic  has  a more  favorable  growth 
rate,  n = 20.1. 
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M ARH  e «s 


t rated  hi  Kig.  11.  Voi  tu  ity  contours  are  shown  at  the  same  times. 


LAYER  DEVELOPMENT  FOR  d = 1/6 , A0=  ,02 


Fin.  19  — Late  development  of  tin*  shear  layer.  For  this  calculation  d = 0.107.  A,,  - 0.02, 
with  fit  - 0.004  sec  for  a 7 X 7 mesh  Gaps  in  the  marker  particles,  particularly  at  the 
free  surface  and  rigid  bottom,  are  artificial  and  are  due  to  round  off  in  tile  low  i recision 
marker  routines. 
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